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1.1  OBJECTIVE 

To  study  active  structure  control  technology  as  applied  to 
Large  Space  Structures  (LSS).  In  particular,  to  develop  analytical 
tools  for  analyzing  their  behavior  and  to  develop  a  means  of 
ensuring  and  testing  for  stability.  This  objective  Is  In  support 
of  the  overall  ACOSS  objective  which  is  to  develop  and  understand 
a  generic,  unified,  structural  dynamics  and  control  technology 
base  for  LSS  with  stringent  Li ne-of-S i ght  (LOS)  and  figure  per¬ 
formance  requirements  that  must  be  maintained  in  the  presence  of 
on-board  and  natural  disturbances. 

1.2  SCOPE 

This  effort  will  be  restricted  to  theory  development,  demon¬ 
stration  and  documentation  of  structural  dynamics  and  control  tech¬ 
nology  for  large  precision  spacecraft.  The  contractor  planned 
the  overall  effort,  defined  the  problem,  and  developed  and  evalu¬ 
ated  solutions.  The  development  and  evaluation  of  solutions 
involved  analytical  studies,  modeling  and  simulations.  The  con¬ 
tractor  identified  and  documented  all  the  analyses,  studies,  de¬ 
signs,  models,  simulations ,  tradeoffs,  issues,  software,  requirements, 
specifi cations  and  data  produced  under  the  contract.  Important 
products  of  the  effort  are  the  .  demonstrati  on  and  documentation  of 
new  control  technology. 
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1.3  STATEMENT  OF  WORK 

The  following  is  a  summary  of  the  tasks  and  technical 
requirements  that  comprise  the  Statement  of  Work  (SOW).  The 
numbers  are  identical  to  those  used  in  the  SOW. 


v4.1.  The  contractor  shall  evaluate  the  existing  LSS  math¬ 

ematical  models  that  have  been  developed  by  Charles  Starke  Draper 
Laboratory  (CSOL).  The  contractor  shall  accomplish  the  following 
tasks: 

4.1.1.  Development  and  Integration  of  Active  Control  Efforts. 

4. 1.1.1.  The  contractor  shall  synthesize  a  high  order  mathematical 
model  of  a  Large  Space  System.  The  model  shall  be  amenable  to 
simplification  (for  purposes  of  analytical  tractab i 1 i ty )  by  trun¬ 
cation.- 

4. 1.1.2.  The  contractor  shall  investigate  sampling  phenomena 
arising  from  a  digital  version  of  an  on-board  controller,  using 
a  truncation  version  of  the  model  developed  under  part  4. 1.1.1. 
above . 

4. 1.1. 3  The  contractor  shall  'integrate  the  technical  approaches 
of  4.1.1,  4.1.2,  and  4.1.3. 

4.1.2  Stability  of  Large  Scale  Systems. 

Using  the  high  order  mathematical  model  of  a  Large 
Space  System  developed  under  4.1.1,  the  contractor  shall  for¬ 
mulate  a  suitable  mathematical  framework.  Using  the  decomposition 
principle,  the  contractor  shall  investigate  its  applicability. 

4.1.3  Development  of  Parameter  Space  Technique. 

The  contractor  shall  develop  the  Parameter  Space 
technique  to  provide  a  desired  set  of  rigorous  stability  criteria 
applicable  to  the  high  order  model  of  the  Large  Space  System  de¬ 
veloped  under  4.1.1.  Its  application  to  the  determination  of 
robustness  will  be  investigated. 

a. 2.  Oral  Presentations.  Oral  presentations  shall  be  given  n 

at  such  times  and  places  as  designated  by  the  Contracting  Officer. 
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1.4  SCHEDULE 


The  original  scheduled  duration  of  the  contract  was  the  period 
1  April  1980  through  31  January  1981.  This  was  subsequently  mod¬ 
ified  to  extend  the  period  to  1  April  1981.  The  reason  was  to 
permit  coverage  of  additional  ACOSS  activities  pending  the  start 
of  a  new  ACOSS  contract. 

1.5  ADMINISTRATIVE 

There  were  four  modifications  to  the  original  contract. 
Modification  #1.  This  was  to  alter  the  originally  estimated 
man-hours  to  comply  with  more  accurate  subsequent  estimates. 
Modification  #2.  This  was  to  change  the  company  address. 
Modification  #3.  This  extended  the  duration  of  the  contract 
to  28  February  1981. 

Modification  #4.  This  extended  the  duration  of  the  contract 
to  1  Apri 1  1 981 . 

1.6  CONTRACT  IMPLEMENTATION 

Dr.  Seltzer  and  Dr.  Worley  evaluated  the  Large  Space  Structure 
mathematical  models. 

Or.  Worley  synthesized  the  "high-  order"  mathematical  model 
of  a  Large  Space  Structure. 

Dr.  Seltzer  investigated  sampling  phenomena  under  the  development 
of  extensions  of  the  Parameter  Space  technique. 
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Professor  Siljak  and  Dr.  Seltzer  investigated  stability  aspec 
of  Large  Space  Structures. 

Professor  Asner,  Dr.  Seltzer,  and  Dr.  Worley  developed  exten¬ 
sions  of  the  Parameter  Space  technique. 
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2.0  EVALUATION  OF  OARPA  MODEL 

2.1  PRELUDE 

The  evaluation  of  any  model  requires  an  understanding  of 
the  requirements  that  caused  the  model  to  be  generated  in  the 
first  place.  The  CSDL  model  was  developed  with  the  stated 
objective  to  assess  the  performance,  sensitivity,  and  hardware 
requirements  of  the  various  active  structural  control  methods 
which  are  being  developed.  In  addition,  the  design  of  the  model 
was  driven  by  the  desire  to  incorporate  certain  attributes  into 
the  overall  system  characteri s ti cs ,  such  as,  (1)  a  structural 
design  based  on  realistic  size  and  weight,  (2)  a  simple  unclas¬ 
sified  optical  system  with  associated  performance  measures  and 
tolerances,  and  (3)  a  set  of  disturbances  typical  of  equipment 
vibration  and  attitude  control.  The  resulting  model  was  pre¬ 
sented  in  Report  RADC-TR-80-377,  "Active  Control  of  Space  Structures 
Interim  Report"  dated  September  1980  by  The  Charles  Stark  Draper 
Laboratory,  Inc.  The  purpose  of  this  section  is  to  review  in 
a  cursory  manner  this  model  to  determine  the  degree  to  which  the 
stated  objectives  are  reached. 

2.1.1  KINETIC  ENERGY 

The  develooment  of  a  structural  model,  or  any  dynamic  model 
for  that  matter,  can  begin  with  the  representation  of  the  kinetic 
and  potential  energy  of  the  system  under  question.  The  general 
form  of  the  kinetic  energy  is  presented  and  then  simplified  under 


the  assumption  that  the  kinetic  enerqy  is  not  a  function  of  the 
generalized  coordinates,  but  rather  only  a  function  of  the  rate 
of  change  of  the  generalized  coordinates. 

T  *  4  y  R  •  R  dm 

R  =  POSITION  OF  ANY  POINT 

T  =  KINETIC  ENERGY 

V  =  POTENTIAL  ENERGY 

d  f  3T  \  3T  3v 

'  3q7  +  *  Qi  ♦  QAi 

ql  =  GENERALIZED  COORDINATE 
Q  =  GENERALIZED  FORCE ( S ) 

QM  2  CONSTRAINT  FORCE(S) 

IF  T  *  T(qi  ) 

THEN 

Tt  3  01  +  QAi 

n  q  y 

2.1.2  PLANAR  EXAMPLE 

We  now  move  from  this  general  form  of  the  kinetic  enerqy 
to  a  planar  representation  of  a  structure  so  that  we  may  be  able 
to  more  completely  understand  the  physics  of  the  situation.  Fo 
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this  purpose,  consider  the  specification  of  the  motion  of  a  structure 
as  depicted.  The  kinetic  energy  and  the  required  oartial  derivatives 
can  be  represented  as  given. 


Pig.  2.1.  Planar  Structure 

LET  r  *  y  +  xe  +  n. 

4>i  3  MOOE  SHAPE 

ni  *  modal  generalized  coordinate 


R  3  r 

J-l  3  j 

WHERE  ♦,  2  t  ,  3  Y 

$2  =  X  n  3  9 


R  3  E  0  n  ;  3R 
j“x  j  j  ’  7? k 


R  dm 


£  4>  &  n  dm 

i*l  J  i  ) 


IF 

3T 

3nk 


j  V* 


M)  "j 


dm  *  M^  6^^ 


SIMILARLY  WITH  POTENTIAL  ENERGY 


3  v 
3lj 


m3  ^ 


u>3  =  jch  EIGEN  VALUE 


THUS 


M4  n  +  M  d>i  n  -  Q.  +  Q* 
3  3  3  J  j  3 


OR  ADOING  "MODAL"  DAMPING 


ALSO  Q  =  F«j(T)  -  F5j 


a  =  SLOPE  AT  F 
3 
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2.1.3  GOVERNING  DIFFERENTIAL  EQUATIONS 

We  now  must  call  on  orthogonality  of  the  "modes"  to 
simplify  the  equations.  Nonorthogonal  modes  can  be  used  and 
are  used  in  some  situations,  but  there  is  increased  coupling 
in  the  equations  with  corres pondi ng  increased  numerical  dif¬ 
ficulties.  The  differential  equations  of  motion  then  follow 
as  given.  It  is  at  this  ooint  that  the  concept  of  equivalent 
viscous  structural  damping  is  introducted  into  the  equations. 

This  representation  of  energy  loss  in  structures  Is  supported 
by  experimental  data  when  the  damping  is  small  (less  then  Zt) , 
but  can  be  seriously  questioned  when  used  for  large  space 
structures . 

2.1.4  DIFFERENTIAL  EQUATIONS  WITH  CLOSED-LOOP  CONTROL 

Now  it  is  possible  to  add  the  control  forces  and  moments 
through  the  concept  of  a  generalized  force.  We  have  taken  here 
a  standard  PD  (position  and  deri vati ve-of-pos i ti on)  controller 
where  the  control  force  or  moment  is  taken  to  be  a  linear  function 
of  the  position  and  the  rate  of  the  body.  Of  course  this  inform¬ 
ation  must  be  sensed,  and  these  control  errors  are  really  only 
an  estimate  of  the  actual  position  and  rate  of  the  body  in  question. 

It  is  now  possible  to  see  the  influence  of  the  control 
force  and/or  moment  on  the  damping  and  stiffness  of  any  given 
structural  generalized  coordinate. 
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2,1.5  SIMPLE  MODEL 

It  Is  now  beneficial  to  simplify  the  model  even  further 
and  consider  a  two-mass  model  as  given  -in  the  illustration.  The 
equations  of  motion  are  simple  to  derive  and  therefore  are  not 
presented.  The  resulting  characteristic  equation  is  presented 
in  a  normalized  form.  The  significant  variables  remaining  in 
the  characteristic  equation  are  the  ratio  of  the  inertias  of  the 
two  bodies,  r,  and  the  ratio  of  the  coupling  of  information 


between  the  two  bodies  in  the  error  signal,  y.  We  now  will 
illustrate  the  effect  of  these  two  parameters  through  the 
use  of  root-locus.  The  numerical  results  are  contained  in 
graphs  found  at  the  end  of  this  section. 


DEFINE 

r  2  I  x  /  I  z 

Y  2  RATIO  OF  BODY  2  TO  BODY  1  INFORMATION 
CONTAINED  IN  ERROR  SIGNAL 

S  3  S/Un 


CHARACTERISTIC  EQUATION 


cigure  2.3  presents  the  root  locus  for  the  case  when  the 
two  bodies  in  our  model  have  equal  moments  of  Inertia.  (Read 
here  that  the  modal  gains  at  each  of  the  two  ooints  on  the 
structure  have  equal  magnitudes  and  opposite  signs.)  The  effect 
of  increasing  the  coupling  between  the  two  bodies  through  the 
control  loop  (increasing)  is  to  decrease  the  damping  in  the 
lightly  damped  mode  and  to  eventually  drive  it  unstable. 

Figure  2.4  is  a  root  locus  for  the  case  when  the  control 
actuator  is  on  a  body  that  hasten  times  the  inertia  of  the 
other  body.  Here  it  can  be  seen  that  the  effect  of  the  coupling 
in  the  control  signal  is  to  drive  the  system  more  quickly  unstable 
The  modal  deflection  at  the  second  point  is  the  culprit,  being 
inversely  proportional  to  the  moment  of  inertia  ratio. 

2.1.6.  OBSERVATIONS 

This  short  journey  into  the  world  of  flexible  body  analysi 
and  control  was  for  the  purpose  of  developing  an  understanding 
of  the  manner  in  which  physical  systems  work  together  and  how 
control  systems  play  in  the  game.  We  believe  that  it  is  very  inroort- 
ant  for  one  to  understand  the  system  on  which  he  is  working 
and  not  just  take  the  most  sophisticated  computer  program  that 
can  be  found  and  start  turning  the  crank.  There  are  a  lot  of 
pitfalls  that  can  be  avoided  if  one  takes  a  lantern  along  with  him 
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2.2  SUMMARY  EVALUATION  OF  CSDL  -  DEVELOPED  MODEL 


This  model  appears  to  have  a  good  representati on  of  the 
optical  bench.  The  complex  paths  that  the  light  rays  follow 
and  the  lengths  of  these  paths  are  sure  to  cause  problems 
if  one  uses  the  focal  plane  error  signal  in  the  control  Iood. 
However,  the  subsystem  structure  is  poorly,  if  not  unsatisfac¬ 
torily,  represented.  In  addition, the  design  of  the  isolator  is 
in  some  question.  The  design  appears  to  be  based  on  the  desire 
to  isolate  longitudinal  vibrations,  while  lateral  and  bending 
vibrations  would  have  a  more  harmful  effect  on  the  image. 

The  following  equations  illustrate  one  type  of  simple 
analysis  that  can  and  should  be  used  to  determine  the  efficacy 
of  any  model.  Here  simple  mathematical  expressions  are  utilized 
to  yield  the  result  that  the  CSDL  Instrument  Section  has  zero 
dimension  in  the  z  direction.  This  is  a  serious  flaw  in  our 
opinion  in  that  this  is  the  section  in  which  the  control  act¬ 
uators  are  present  and  also  the  system  rate  gyros.  We  have 
already  illustrated  that  flexibility  between  these  two  devices 
is  a  serious  problem  for  control  system  designers. 
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3.0  CONTROL  DYNAMICS  MODEL 
3.1  ANALYSIS 

The  purpose  of  this  section  is  to  present  an  alternate 
approach  to  the  model  development.  We  presented  a  model  to 
DARPA  May  8,  1980  in  which  we  felt  were  all  the  ingredients 
to  properly  cause  a  control  system  designer  to  work  for  his 
lunch.  This  model  was  derived  using  the  transfer  matrix 
technique  based  on  an  outline  presented  in  "Analytical  Methods 
in  Vibrations"  by  Leonard  Meirovitch.  Press  of  time  and  space 
will  not  allow  a  complete  description  of  the  development  pro¬ 
cedure,  but  the  complete  modal  model  was  developed  using  a 
Hewlett  Packard  desktop  comDuter,  Model  9835A. 

A  picture  of  the  model  is  presented  on  Fig.  3.1.  Here 
we  have  attempted  to  include  all  the  features  that  would  cause 
the  control  system  designer  to  have  troubles  in  coming  up  with 
a  stable  and  wel 1  - performi ng  system.  The  mode  shapes  arising 
out  of  this  model  are  presented  in  the  fallowing  pages. 

The  first  bending  mode  (Fia.  3.2.)  is  being  driven  by  the 
stiffness  of  the  interface  and  the  transverse  moment  of  inertia 
of  the  optical  bench  and  the  equipment  section.  For  this  model, 
the  second  mode  (Fig.  3.3.)  is  primarily  a  torsional  mode  with 
the  natural  frequency  being  driven  by  the  torsional  stiffness  of 
the  interface  and  the  roll  moment  of  inertia  of  the  two  sections. 
The  third  mode  (Fig.  3.4.)  is  a  combined  bending  and  torsion  mode 
of  the  optical  bench.  The  mode  shape  and  frequency  of  this  mode  is 
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determined  by  the  properties  of  the  optical  bench  alone  .  The 
fourth  mode  (Fig.  3.5.)  is  a  pure  bending  (first  mode)  of  the 
optical  bench  and  as  such  is  completely  determined  by  the  physical 
properties  of  that  section.  The  fifth  mode  (Fig.  3.6.)  is  a 
complex  mode  involving  first  bending  of  the  optical  bench  with 
flexing  of  the  interface  structure.  It  would  be  difficult  to 
characterize  data  to  change  this  mode,  but  one  would  expect 
stiffness  of  the  interface  to  have  some  effect.  The  sixth 
mode  (Fig.  3.7.)  represents  a  second  bending  of  the  optical 
bench  and  as  such  is  a  brother  of  mo^e  #4.  The  seventh  mode 
(Fig.  3.8.)  is  pure  torsion  of  the  equipment  section.  It  also 
represents  graphically  the  effect  of  the  isolator  section  in  that 
it  is  completely  uncouoled  from  the  optical  bench.  The  eighth 
mode  (Fig.  3. 9. ), 1  ike  mode  #7,  is  completely  uncoupled  by  the 
isolator  section  and  represents  pure  second  mode  torsion  of  the 
optical  bench. 

3.2  CONCLUSIONS 

We  have  presented  the  results  of  a  simple  model  of  an  example 
large  space  structure.  The  results  are  in  the  form  of  mode 
shapes  and  natural  frequencies  as  presented.  In  addition,  because 
of  the  simplicity  of  the  model  one  is  able  by  inspection  to 
ascertain  the  significant  configurational  drivers  for  each  of  the 
modes  and  frequencies  and  thereby  be  in  the  position  to  change 
them  in  an  appropriate  manner,  thereby  providing  "management 
flexibility."  It  is  through  this  process  of  developing  and  main- 
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talnlng  an  understanding  of  not  only  the  control  process  but  also 
the  open  loop  plant  charac terl s ti cs  that  we  believe  is  the  most 
effective  way  to  design  control  systems  for  vehicles.  Therefore, 
it  is  also  believed  that  if  one  is  attempting  to  design  a 
structure  that  would  exercise  control  theorists,  then  one  should 
pursue  a  similar  course.  We  did  not  always  have  these  qlant 
computer  programs  with  us  and  large  SDace  structures  were  flown 
very  successfully  in  the  past  (for  example,  the  SATURN  IB  and 
V  launch  vehicles  and  the  SKYLAB  manned  satellite).  Finally, 
of  course,  the  program  costs  and  time  are  minimal. 


4.0  VALIDATION  OF  REDUCED-ORDER  MODELS  FOR  CONTROL  SYSTEM  DESIGN 
The  concept  of  suboptimality  is  applied  to  testing 
validity  of  reduced-order  models  in  design  of  feedback,  schemes 
for  large-scale  systems.  Aggregation  and  singular  value  decom¬ 
position  as  model  reduction  techniques,  are  interpreted  in 
the  expansion-contraction  framework,  which  is  suitable  for 
evaluation  of  suboptimality  of  closed-loop  systems  resulting 
from  reduced  order  designs.  The  proposed  validation  procedure 
is  applied  to  a  control  design  of  a  large  space  structure. 


4.  1  INTRODUCTION 


In  system  theory  and  applications,  "model  reduction"  stands  for  a 
variety  of  more  or  less  disparate  concepts  and  techniques  which  have  a 

common  goal:  to  reduce  the  size  of  the  mathematical  model  of  a  large 

\ 

system  in  order  to  sinplify  the  design  of  control  and  estimation  schemes. 
IXie  to  the  nature  of  the  mathematical  models  involved  in  describing  the 
motion  of  Large  Space  Structures  (LSS),  a  large  number  of  recent  results 
in  model  reduction  have  been  obtained  in  this  field,  with  enphasis  on  the 
size  and  complexity  of  the  controller  to  be  implemented  onboard  the  space¬ 
craft  (see,  for  exanple,  fl]  and  [2]).  Most  of  the  proposed  model  reduc¬ 
tion  procedures  has  a  common  plan:  First,  the  elements  of  a  mathematical 
model  are  separated  into  two  parts  according  to  the  significance  of  their 
contribution  to  the  desired  performance  characteristic  of  the  system  (sta¬ 
bility,  optimality,  frequency  response,  and  the  like).  Second,  the  most 
significant  part  of  the  model  is  used  to  design  the  control  scheme  for  the 
overall  system.  Finally,  the  scheme  is  implemented  and  the  system  perfor¬ 
mance  is  tested  usually  by  an  extensive  simulation  analysis.  The  enphasis 
on  each  phase  of  the  model  reduction  plan  varies  from  concept  to  concept, 
but  the  first  phase  involving  the  decorrposition  of  the  model  has  been  the 
most  considered  area  in  the  research  on  model  reductions. 

Almost  exclusively,  all  the  model  reduction  techniques  are  concerned 
with  the  open-loop  behavior  of  the  system  and  are,  therefore,  subject  to 
the  question  of  whether  the  reduced-order  model  would  still  bo  a  go xi  n<- 
presentation  of  the  overall  system  aft r,r  the  feedback  loop  is  closed.  The 
answer  to  this  question  would  nost  probably  be  positive  if  a  reduo 'd-order 
model  were  indeed  a  good  approximation  of  the  original  system.  Yet,  it  is 
difficult  to  say  how  good  the  approximation  is  without  a  ccrrparison  of  the 


resulting  closed-loop  system  with  the  one  that  would  be  obtained  under  ideal 
conditions  when  a  full  controller  is  designed  iKing  the  entire  model  of  the 
system. 

In  order  to  resolve  the  validation  problem  of  model- reduct ion  schemes, 
we  note  that  there  are  two  distinct  simplifications  attempted  by  the  nethods: 
To  obtain  a  simple  control  law  which  is  suitable  for  implementation  and  to 
reduce  the'  ccrputations  involved  in  the  control  design.  Since  the  off-line 
computations  needed  to  determine  the  best  control  law  for  dynamic  systems  of 
considerable  size  may  well  be  within  the  capabilities  of  today's  computers, 
the  dimensionality  problem  is  not  as  severe  as  the  compatibility  problem  of 
the  control  structure  with  the  design  constraints.  For  this  reason,  an  ap¬ 
propriate  approach  to  the  model  reduction  problem  is  to  concentrate  on  con¬ 
struction  of  suitable  control  laws  and,  subsequently ,  compare  the  resulting 
closed- loop  system  with  the  optimal  off-line  model  used  as  reference.  This 
approach  to  model  reduction  is  in  the  spirit  of  subopt imal  design  schemes 
developed  in  the  context  of  large-scalp  systems  [3,  4],  and  the  major  objec¬ 
tive  of  this  paper  is  to  adapt  these  suboptimality  schemes  for  use  in  model- 
reduction  problems.  This  reformulation,  in  turn,  opens  up  a  real  possibility 
for  future  use  of  the  model- reduct ion  approach  in  the  design  of  modem  decen¬ 
tralized  control  schemes  for  large  space  structures  composed  of  interconnected 
systems. 

The  organization  of  this  section  is  as  follows: 

In  the  Section  4.2  ,  we  briefly  summarize  the  concept  of  suboptimality 
which  provides  a  basis  for  the  subsequent  development .  In  Section  4.3,  we  out¬ 
line  the  inclusion  principle,  which  const  itutes  an  appnpriatc  mathemat  ical 
framework  for  perfect  model  reduction.  In  this  section,  we  also  study  the 
role  of  inclusion  principle  in  the  context  of  suboptimal  Control,  and  show 
that  when  perfect  model  reduction  is  possible,  then  the  suboptimal  control 
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obtained  from  the  reduced-order  model  is  actually  optimal  for  the  original 
overall  system.  In  Section  4.4,  we  elaborate  on  various  approximate  model 
reduction  schemes.  From  the  point  of  view  of  suboptimalitv,  these  schemes 
appear  as  deviations  from  perfect  model  reductions  obtained  via  the  inclu¬ 
sion  principle.  In  Section  4.5,  a  model  reduction  technique  based  upon  the 
balanced  representations,  is  applied  to  a  model  of  a  large  space  structure. 
Finally,  in  conclusion,  we  discuss  possible  extensions  of  the  suboptimal 
approach  to  more  elaborate  problems  of  LSS. 
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4.2  SUBOPTIMAL  CONTROL 


Let  us  consider  a  linear  system  S  described  as 


P 


S:  x  -  Ax  +  Bu  ,  x(0)  »  xQ  (2.1) 

wnere  x(t)  6  ^  is  the  state  and  u(t)  6  fiP  is  the  input  of  S  ,  and  the 
matrices  A  and  B  are  constant  and  of  appropriate  dimensions.  With  S  we 
associate  a  quadratic  performance  index 

OD 

J  *  |  (xTQx  +  uTRu)dt  ,  (2.2) 

0 


where  Q  is  a  symmetric  nonnegative  definite  matrix  and  R  is  a  symmetric 
positive  definite  matrix,  both  matrices  being  constant  and  of  appropriate 
dimensions.  The  system  5  and  the  cost  functional  J  form  a  standard  opti¬ 
mization  problem  {S,  J}  ,  the  solution  of  which  is  given  as 

u°=-K°x  ,  (2.3) 

where  the  gain  matrix  K°  is 

K°  -  R-1B*P  ,  (2.4) 


and  P  is  the  solution  of  the  Riccati  equation 

A^P  +  PA  -  PBR-1BP  +  Q  =  0  .  (2.5) 

Application  of  the  control  u°  to  S  results  in  the  optimal  closed-loop 
system 

S° :  x  =  (A-BK°)x  ,  (2.6) 

and  yields  the  optimal  cost  J°(Xq)  =  K°)  , 

J^Xq)  -  xJ’Xq  .  (2.7) 
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New,  suppose  instead  of  the  optimal  control  u  ,  we  apply  another  con- 

u  -  -to  ,  (2.8) 

where  the  gain  matrix  K  is  obtained  by  some  rational  procedure  which  is 
different  from  the  above  optimal  solution  of  {S,  J}  .  Then,  the  control  u 

of  (2.8)  produces  a  cost  J  (xQ)  *  J(j^,  K)  ,  which  is,  in  general,  larger 

o  • 

than  the  optimal  cost  J  (Xq)  .  In  fact,  the  cost  J  (Xq)  may  very  'veil  be 
infinite  if  the  control  in  (2.8)  is  not  stabilizing.  We  now  state  the  follow¬ 
ing  [  4  ] : 

(2.9)  DEFINITION.  The  control  lots  (2.8)  ia  said  to  be  suboptimal  with 
degree  y  for  S  with  respect  to  S°  if  there  exists  a  positive  number  y 
such  that 

J*(xQ)  <  y'ViXjj)  ,  (2.10) 

for  all  Xq  . 

lb  derive  conditions  for  suboptimality  and  conpute  the  suboptimality 
index  y  ,  we  note  that  if  the  control  (2.8)  is  stabilizing  for  S  ,  then 
the  suboptimal  cost  is  given  by  • 

J#(x0)  -  xjllx0  ,  (2.11) 

wfaere  H  is  the  unique,  symmetric,  and  positive  definite  solution  of  the 
Liapunov  equation 

(A-DK)T!I  +  II(A-Hi)  +  Q  +  KRK  =  0  .  (2.12) 

This  leads  to  [ 5  ] : 

(2.13)  THEOREM.  If  the  control  law  (2.8)  is  stabiliziiig  for  S  ,  then  it  is 
suboptimal  for  S  with  degree 
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where  denotes  the  largest  eigenvalue  of  the  indicated  matrix. 


An  alternative  suboptimality  criterion,  which  does  not  require  a  test 
for  stability,  is  provided  by  the  following  [4  ]: 

(2.15)  THEOREM.  The  control  law  (2.8)  i8  suboptimal  for  S  with  degree 
U  if  the  matrix 

F(u)  -  (K°-K)TR(K°-K)  -  (l-y)(Q  +  KTRK)  (2.16) 

ia  nonnegative  definite. 

We  note  that  to  confute  the  suboptimality  index  u  using  Theorem  (2.13) 
or  (2.15),  we  have  to  conpute  the  solution  of  the  problem  {S,  J)  .  Although 
we  know  the  optimal  control  (2.3),  we  may  prefer  to  apply  the  suboptimal  con¬ 
trol  (2.8).  One  irrportant  reason  might  be  that  the  optimal  control  is  diffi¬ 
cult  to  irrplement  in  the  actual  design.  For  exanple,  we  may  choose  the  output 
feedback 

u  *  — PV  ,  (2.17) 

where 

y  =  Cx  (2.18) 

is  the  only  available  measurement  with  the  matrix  C  having  a  full  row  rank. 
In  this  case,  a  reasonable  choice  for  K  would  be 

K  -  K°CT(CCT)-1C  ,  (2.19) 

which  is  the  projection  of  K°  on  the  row  space  of  the  output  matrix  C  . 

We  do  not  gain  anything  computationally ,  but  only  satisfy  the  design  con¬ 
straints  on  the  control  structure.  Another  way  to  generate  the  gain  matrix 
K  in  (2.8),  which  satisfies  the  information  structure  constraint,  is  to  re- 


place  by  zeros  those  elements  of  K°  that  correspond  to  states  which  are 
not  available  for  feedback.  Following  [  4  J ,  both  types  of  control  are  re¬ 
ferred  to  as  degenerate  control. 

Alternatively,  the  control  (2.8)  can  be  generated  frcm  the  solution  of 
a  smaller  problem  {5,  3}  where  the  system 

3:  x  *  Ax  +  Bu  ,  x(0)  ■  SL  (2.20) 


represents  a  reduced  order  model  of  S  ,  that  is,  x(t)  e  and  n  _<  n 


The  cost 


J  =  [  (xTQx  +  uRu)dt  , 


(2.21) 


is  obtained  from  the  cost  J  in  a  nanner  conpatible  with  the  reduction  of 
S  to  3  .  Usually,  S  can  be  chosen  as  a  projection  of  S  onto  sane  sub¬ 
space  of  the  state  space  [5,6]  defined  by  a  singular  transformation 


x  *  Ux  . 


(2.22) 


Now,  if  the  optimal  solution  of  the  problem  (3,  3}  is  obtained  as 
u  ■  -R°x  , 


(2.23) 


then  the  control 


u  =  -K°Ux 


(2.24) 


can  be  a  reasonable  choice  for  the  original  system  S  .  Suboptimality  of  the 
control  (2.24)  can  be  tested  using  Theorem  (2.13)  or  (2.15)  with  K  *  R°U  . 

Another  advantage  of  obtaining  the  suboptimal  control  by  the  solution 
of  a  smaller  optimization  problem  {S,  3)  ,  is  a  possibility  to  oonpare  the 
suboptimal  cost  J*(Xq)  *  J(xQ,  K°U)  with  the  optimal  co6t  3°(x0)  *  3(Xq,  R°) 
of  pair  {3,  3)  .  Of  course,  should  be  chosen  as  the  projection  of 
Xq  ,  that  is,  Xq  ■  Uxq  .  If  this  corrparison  can  be  worked  out,  then  we  can 


avoid  solving  the  original  problem  {S,  J}  ,  thus  providing  for  both  conpu- 
tational  and  inplementational  sinplifications.  In  the  context  of  suboptimal¬ 
ity  [  4  ,  5  ] ,  toldjig  the  optimal  small  system  S°  as  reference  in  the  defi¬ 
nition  of  suboptimality ,  is  the  dual  of  taking  the  original  system  S°  as  the 
reference  system,  in  the  sense  that  in  the  latter  case  the  control  is  considered 
to  be  perturbed  from  the  optimal  value,  while  in  the  former  case  the  system  is 
perturbed  by  the  residuals  which  were  eliminated  by  the  projection. 
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4.3  THE  INCLUSION  PRINCIPLE 


In  using  a  reduced  order  model  to  build  a  controller,  it  is  important  to 
know  if  the  behavior  of  the  original  system  can  be  predicted  by  the  rodel. 

The  inclusion  principle  of  [5,6  ]  provides  conditions  under  which  the  re¬ 
duced  order  model  represents  perfectly  the  motion  of  the  original  system.  In 
this  section,  we  briefly  sunrarize  the  principle  and  discuss  its  triplications, 
let  us  reconsider  the  two  systems  S  and  3  ,  which  are  described  now  as 

S:  x  *  Ax  +  Bu  ,  x(0)  =  Xq 

y  •  Cx  (3.1) 

and 

3:  x  =  Ax  +  Bu  ,  x(0)  - 

y  •  Cx  (3.2) 

where  the  outputs  y(t)  €  R*  and  y(t)  €  R*  are  added  to  the  previous  des¬ 
criptions  (2.1)  and  (2.20).  The  solutions  of  (3.1)  and  (3.2)  are  denoted  by  . 
x(t,  xQ;  u)  and  x(t,  5^;  u)  ,  respectively.  Ife  have: 

(3.3)  EEFTNOTON.  The  system  S  of  (3.1)  is  said  to  include  the  system  3 
of  (3.2)  if  there  exist  matrices  U  and  .V  of  order  n  x  n  and  n  *  5  such 
that  UV  =  I_  and 


x(t,  XqJ  u)  -  Ux(t ,  VxQ;  u) 

(3.4) 

y(t)  *  y(t) 

(3.5) 

for  all  t,  Xq  ,  and  u(t)  . 

Definition  (3.3)  inplies  that  if  S  includes  3  ,  then  the  system  3 
contains  all  the  information  about  the  system  S  .  If  S  includes  3  ,  then 
S  is  said  to  be  an  expansion  of  3  .  and  3  is  a  contraction  of  S  .  In 


(3.3),  I-  denotes  the  n  *  n  identity  matrix.  Obviously,  an  interesting 
question  is:  Under  what  conditions  S  includes  2  ?  This  is  provided  by 
the  following: 

(3.6)  THEOREM.  S  includes  2  if  and.  only  if  there  exist  matrices  U  and 

V  such  that  UV  *  I-  and 
n 

UAiV  ■  A1,  UA1_1B  -  Si_1B,  CA1-1V  *  CA1”1,  CAi-1B  *  CAi_1B 

i  -  1,  2,  ...  ,  n  .  (3.7) 


Let  us  consider  some  special  cases.  First,  suppose  there  exists  an 
n  x  n  matrix  U  with  full  row  rank  such  that 

UA-SU,  UB  *  B  ,  C-&J  .  (3.8) 

Then,  it  is  easy  to  show  that  the  conditions  in  (3.7)  are  satisfied  for  any 
right  inverse  V  of  U  ,  that  is,  for  any  V  such  that  UV  *  I-  .  Therefore, 
if  the  conditions  (3.8)  are  satisfied  for  some  U  ,  then  2  is  a  contraction 

t 

of  S  and,  hence,  (3.4)  and  (3.5)  hold.  But,  in  this  special  case,  there  is 
a  stronger  relation  between  the  solutions  of  S  and  2  ,  namely, 

x(t,  UXq;  u)  -  Ux(t,  xQ;  u)  ,  (3.9) 

for  all  t,  Xq  ,  and  u(t)  .  In  other  words,  whatever  the  initial  state  xQ  , 
a  projection  of  the  solution  of  S  can  be  recovered  from  a  solution  of  2  . 

In  this  case,  2  is  an  aggregation  of  S  in  the  sense  of  Aoki  [7  ]. 

As  another  special  case,  suppose  there  exists  an  n  *  n  matrix  V  with 
full  oolunn  rank  such  that 

AV  «  VX  ,  B  -  VO  ,  CV  =  C  ’.  (3.10) 

Then,  for  any  left  inverse  U  of  V  satisfying  UV  =  1-  ,  conditions  (3.7) 
hold,  and  2  is  a  contraction  of  S  .  In  this  case,  too,  a  stronger  relation 
holds, 
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Vx(t,  xQ;  u)  -  x(t,  V5cQ;  u) 


(3.11) 


which  inplies  that  solutions  of  S  starting  in  a  subspace  of  the  state  space 
of  S  can  be  recovered  from  the  solutions  of  3  .  In  this  case  we  say  that 
3  is  a  restriction  of  S  . 

Let  us  now  broaden  the  inclusion  principle  to  consider  the  optimization 
problems  via  model  reduction.  In  the  optimal  control  problem  (S,  J)  ,  we 
let 

Q  =  CTC  (3.12) 


for  -some  matrix  C  and  rewrite  J  as 


|  (yTy  +  uTRu)dt 


(3.13) 


where  y(t)  may  be  taken  as  the  output  of  S  described  by  (3.1).  Let  us 
assure  that  the  system  5  of  (3.1)  is  aggregable  to  a  system  3  described 
by  (3.2).  With  3  we  associate  the  cost 


j  (yTy  +  uTRu)dt  . 


(3.14) 


We  choose  a  control 
u  -  -JS  , 


(3.15) 


which  is  not  necessarily  optimal  for  (3,  J}  ,  and  apply  to  3  .  Lhder  the 
assumption  that  it  is  stabilizing,  the  oontrol  (3.15)  yields  the  cost  given 
by 

3(x0,  K)  =  xjnx0  (3.16) 

where  H  is  the  unique,  synmetric,  and  positive  definite  solution  of  the  Lia¬ 
punov  matrix  equation 

rp 

(3.17) 
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(A-BR)Th  +  R(S-ER)  +  CTC  +  RtkR  -  0  , 


Multiplying  (3.17)  by  U  from  the  left  and  by  D  from  the  right,  and  using 


the  aggregation  conditions  (3.8),  we  get 


where 


(A-BK)TH  +  H(A-BK)  +  CTC  +  KTRK  -  0  , 


k  *  Ru  ,  h  -  irfiu  . 


(3.18) 


(3.19) 


Comparing  (3.18)  with  (2.1 2),  and  keeping  in  mind  that  Q  *  C1C  ,  we  con¬ 


clude  that 


RU)  -  xJ^HUXq  -  3(1^,  R) 


(3.20) 


In  particular,  if  R  *  R°  ,  the  optimal  feedback  matrix  for  the  problem  {5,  3}  , 


then  (3,20)  gives 


J(xQ,  R°U)  -  3(0^,  R°)  -  3°(Ux0)  . 


(3.21) 


Since  3°(UXq)  is  the  global  minimum  of  3  ,  (3.21)  further  inplies  that 
R°U-K°,  (3 


J(Xq,  R°U)  -  J°(Xq)  -  3°(Ux0)  , 


(3.22) 


(3.23) 


which  establish  the  relation  between  the  optimal  controls  and  costs  of  the  two 
optimization  problems  (S,  J}  and  (3,  3)  when  3  is  an  aggregation  of  S  . 
This  result  could  have  been  achieved  directly  by  comparing  the  Riccati  equa¬ 
tions  corresponding  to  the  two  problems.  We  preferred  the  above  derivation, 
because  we  can  reproduce  it  in  the  suboptimality  context  outlined  in  the  next 


section. 


Similarly,  if  3  is  a  restriction  of  S  ,  we  can  show  that 
J(V^,  K)  =  3(3^,  KV)  , 


(3.24) 


provided  that  costs  are  finite.  In  particular, 


w 


4.4  SUBOPT IMAL ITY  OF  REDUCED  ORDER  MODELS 

Che  of  the  difficulties  in  applying  the  inclusion  principle  to  model 
reduction  is  the  fact  that  it  may  be  overly  restrictive  and  conditions  for 
a  given  system  S  to  have  a  contraction  3  do  not  hold.  A  natural  way  to 
resolve  this  problem  is  to  introduce  an  approximate  contraction  procedure 
and  evaluate  the  suboptimality  of  the  end  result. 

It  is  clear  from  the  aggregation  conditions  (3.8)  that  the  system  S 
is  aggregable  to  a  smaller  system  3  if  and  only  if  it  is  unobservable. 

Thus,  if  the  problem  (S,  J}  is  such  that  all  the  state  variables  are 
observed  in  the  performance  index  J  ,  then  it  cannot  be  reduced  to  a 
smaller  problem.  In  order  to  produce  an  approximate  aggregation,  we  split 
the  system  matrix  A  into  two  matrices  as 

A  =  Ajj  +  Aj,  .  (4.1) 

By  choosing  the  "surplus  matrix"  Ag  appropriately,  the  "new"  matrix  A^A-Ag 
can  be  made  such  that  the  system 

SN:  x  -  AjjX  +  Bu 

y  =  Cx  (4.2) 

is  aggregable  to  a  smaller  system 
3:  4  =  Ax  +  Su 

y  =  CR  .  (4.3) 

The  aggregation  conditions  (3.8)  require  that 

U( A  -  Ag)  =  AU  ,  C  =  £U  ('K',) 

for  some  U,  5  ,  and  5  .  The  choice  of  U  is  arbitrary  provided  the  second 
equation  in  (4.4)  is  satisfied  by  soitp  C  .  An  obvious  choice  is  U  =  C  , 

C  ■  I  ,  which  produces  3  of  smallest  dimension.  This  may,  however,  require  * 
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the  elements  of  A  to  be  large  conpared  to  those  of  A  ,  which  is  not  de- 

O 

sirable  since  we  want  to  make  3  "close"  to  S  .  Alternatively,  making  uso 
of  the  concept  of  chain  aggregation  [  8  ]  ,  we  can  choose  U  -  [CT  l£]T, 

v 

C-  tl  0]  ,  \»*iere  is  such  that  U  has  full  row  rank.  In  this  case,  a 

suitable  choice  of  Ur  can  make  it  possible  to  satisfy  (4.4)  with  A  suffi- 

S 


ciently  small. 

With  U  fixed,  from  (4.4)  we  get 

**  UA  -  Au  ,  (4.5) 

which  suggests  that  A  should  be  chosen  as 

A  -  UA  U^UUV1  (4.6) 

in  order  to  minimize  |  jUAg|  |  .  With  A  as  in  (4.6),  (4.5)  reduces  to 

UAg  •  UA  [I  -  UT(UUT)'1U]  ,  (4.7) 

the  minimum-norm  solution  of  which  is 

Ag  =■  UT(UUT)"1UA  [I-UT(UUT)"1U]  .  (4.8) 


In  surmnary,  we  showed  that  by  subtracting  a  matrix  A-.  from  the  nominal 
matrix  A  ,  we  can  produce  an  approximate  aggregate  J  of  S  for  some  U 
and  C  .  The  matrix  A  of  3  can  be  determined  from  (4.6)  and  B  *  UB  . 

We  can  then  proceed  to  solve  the  reduced  order  problem  (3,  3}  and  produce 
the  low-order  gain  matrix  1?  .  Then,  the  optimization  problem  (S  ,  J}  can 
be  solved  to  serve  as  a  reference  for  suboptimality  of  the  system  S  driven 

by  the  control  u  =  R°Ux  of  (2.24),  which  is  obtained  from  the  low-order 
design  R°  .  Subcptimality  conditions  are  provided  by  Theorem  (2.13)  or 
(2.15). 
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to  alternative  to  the  approximate  aggregation  considered  above  is  to 
aggregate  the  weighting  matrix  C^C  of  the  state  of  S  in  the  performance 
index  J  ,  as  proposed  by  Aoki  [7].  This  is  done  by  choosing  an  aggregation 
matrix  U  with  full  row  rank,  such  that 

UA  »  AU  C4*9) 

for  some  matrix  5  ,  and  conputing  B  as 

g  =■  UB  .  (4.10) 

Conditions  (4.9)  and  (4.10)  are  part  of  the  aggregation  conditions  in  (3.8). 

Itxvever,  since  we  assumed  that  S  is  not  perfectly  aggregable  to  a  smaller 
system,  there  is  no  C  that  satisfies  the  last  aggregation  condition,  namely, 
C  »  Cu  .  Here,  we  deviate  from  perfect  aggregation,  and  choose  C  such  that 
|  |C-Cu|  |  is  minimum,  i.e. , 

t  -  CuVf1  .  (4.11) 

The  matrices  {5,  B,  C}  of  (4.9)-(4.11)  constitute  an  approximate  aggregate 
system  S  of  (4.4),  with  which  we  associate  the  cost  7  in  (3.15). 

Nav,  consider  the  Riccati  equation 

A^P  +  PA  -  PBR_1b"IP  +  CTC  =  0  (4.12) 

T 

associated  with  the  aggregate  problem  {5,  7}  .  Multiplying  (4.12)  by  U 
from  the  left  and  by  U  from  the  right,  and  using  (4.9  )-(4Jl  )  we  obtain 

at(utPu>  +  (utpu )a  -  (utpu)br-1bt(utpu)  +  ncTcn  =  o  (4.13) 

where  the  matrix 

n  =  UT(UUT)-1U  ,  (4.14) 

is  a  projection  on  R0  .  Conparing  (4.13)  with  the  Riccati  equation  (2.5) 
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T 

of  the  original  problem  {S,  J}  ,  we  observe  that  U  PU  can  be  considered 
as  an  approximate  of  P  provided  CTI  is  close  to  C  .  Thus , 

K  -  R°U  »  R-1BTPU  *  R-1BT(UTPU)  (4.15) 

is  an  approximate  to  the  optimal  feedback  K°  in  (2.4),  the  suboptimality  of 
which  can  be  tested  using  Theorem  (3.13)  or  (3.15). 

Ufe  note  that  condition  (4.  9  )  requires  that  U  consist  of  reciprocal 
eigenvectors  of  A  ,  and  A  retain  those  eigenvalues  of  A  that  correspond 
to  these  eigenvectors.  Thus,  the  aggregation  of  the  performance  index  as 

explained  above  consists  of  choosing  an  eigenspace  of  A  and  modifying  the 
performance  index  to  reflect  the  contribution  of  the  selected  modes  to  the 
cost.  In  this  sense,  this  type  of  approximate  aggregation  is  closely  related 
to  Skelton's  model  reduction  scheme  through  modal  cost  analysis  [9  ,10].  For 
a  given  order  n  of  the  reduced  order  model,  the  n  modes  that  contribute 
most  to  the  cost  J  can  be  determined  by  computing  U  such  that  1 1 C-CII 1 1 
is  minimum  subject  to  (4.9  ),  where  n  is  given  by  (4.14).  However,  whether 
the  nodes  determined  by  this  U  matrix  will  be  the  most  significant  ones  in 
the  closed- loop  system  is  not  cle'ar.  Che  way  of  deciding  on  the  best  approxi¬ 
mate  aggregation  might  be  to  oonpute  the  degree  of  suboptimality  u  for  each 
of  the  n!/(n-n)!R!  possible  choices  of  U  ,  and  to  choose  the  one  with  mini¬ 
mum  u  .  However,  noting  that  the  contribution  of  a  mode  to  the  closed- loop 
cost  depends  on  how  controllable  as  well  as  how  observable  it  is,  a  more  re¬ 
fined  model  reduction  procedure  can  be  developed  in  the  suboptimality  framework 
via  singular-value  deccrposition. 

The  singular-value  decomposition  of  a  system  using  balanced  representation 
has  been  developed  by  Moore  {11]  and  later  studied  in  detail  by  others  (e.g. , 
see  [12,  13]).  In  simple  terms,  a  balanced  system  is  one  which  is  as  control- 
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lable  as  it  is  observable,  where  the  mensure  of  controllability  and  observ¬ 
ability  is  provided  by  the  singular  values  [11]  of  the  respective  grairminns. 
To  illustrate  the  idea,  suppose  S  is  stable  and  let 


and 


Ec  =  |  eAtBBTeATtdt  , 
0 

Iq  =  |  eA  tCTCeAtdt  . 
0 


(4.16) 


(4.17) 


denote  the  controllability  and  observability  granmians  of  S  respectively. 

Now  consider  the  product  2Q2C  .  Assuming  S  is  controllable  and  observable 

(otherwise ,  a  preliminary  procedure  can  be  used  to  reduce  the  system  order  by 

deleting  uncontrollable  and/or  unobservable  part),  2Q  and  2C  are  positive 

definite  matrices,  and  therefore,  all  the  eigenvalues  of  2^-  are  positive. 

2 

Let  these  eigenvalues  be  denoted  by  with  _>  o2  _>  ...  aR  >0  .  The 
positive  square  roots  of  the  eigenvalues  of  2Q2C  are  called  the  singu¬ 
lar-values  of  the  system  S  . 

It  was  shewn  by  Moore  [ll]  that,  there  always  exists  a  coordinate  frame 
in  which 

2C  *  20  "  1  "  ^i»  c2'  »  °n)  *  (4.18) 


The  representation  of  S  in  this  coordinate  frame  is  said  to  be  balanced. 
Now,  suppose  the  system  S  has  a  balanced  representation,  and  let 

2  =  diag  (21#  I2)  ,  (4.19) 

where 

^ 2  =  ^*^6  ***  *  °n ^  * 

22  =  diag  (o.+1,  a.+2 . on)  .  (4.20) 

Let  us  consider  the  corresponding  partition  of  {A,  B,  C)  : 
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A  - 

A11  A12 

,  B  - 

Bl" 

J*21  A22_ 

% 

C  * 

1 — 1 

• 

(4.21) 


If  a.  »  a..,  in  (4.20),  then  we  can  consider  {A,,,  B, ,  C, }  as  the  most 

n  n+i  l  l 

controllable  and  observable  part  of  the  system  S  .  Thus,  the  n-th  order 
system  5  =  (5,  B,  C)  with  3  ■  A^,  B  =  B^,  C  *  ,  can  be  taken  as  a 

good  reduced  order  model  for  S  .  In  fact,  it  is  easy  to  see  from  (4.19) 
that  3  is  a  gpod  approximation  to  S  in  the  sense  that  the  differences 
between  the  controllability  and  observability  gramnians  of  3  and  S  are 
minimum  in  the  norm  square  sense. 

Once  a  reduced  order  model  3  of  the  system  S  is  obtained  through 
the  singular-value  decomposition  of  a  balanced  representation  of  S  ,  a 
subopt imal  control  can  be  designed  for  S  following  the  procedure  of  Sec¬ 
tion  2.  This  approach  to  model  reduction  is  particularly  promising  since 
it  allows  for  generating  systematically  a  sequence  of  reduced  order  models 
consisting  of  smaller  and  smaller  portions  of  the  A,  B,  C  matrices  in 

(4.21).  Computing  the  degree  of  suboptimality  at  each  step  one  can  com¬ 
promise  between  the  degree  of  suboptimality  and  the  order  of  the  reduced 
model. 

Finally  we  note  that  order  reduction  through  singular-value  deconposi- 
tion  can  be  considered  as  a  combination  pf  the  previous  two  approaches, 
where  we  first  subtract  from  the  A  matrix  in  (4.21)  the  surplus  matrix 


(4.22) 

so  that  with  U  =  [I  0]  ,  U(A  — A~)  =  A  U  ,  and  then  choose  5  as  in  (4.11  ), 
n  o 

which  yields  C  =  .  However,  this  procedure  is  not  justified  unless  S  is 


balanced. 
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4.5  APPLICATION 


To  illustrate  the  validation  of  model  reduction  scherres  using  the  con¬ 
cept  of  subopt imality,  we  consider  the  feedback  design  of  a  large  space  struc 
ture  described  in  [14],  The  mathematical  model  is  represented  as  a  set  of 
nine  second  order  equations 

+  d^  +  o^zi  *  bTu  ,  i  »  1,  2,  ...  ,  9  ,  (5.1) 

where  z^t)  are  the  vibration  modes  with  damping  dA  and  frequency  u>.  , 

T 

and  u(t)  is  the  two  dimensional  control  input  such  that  b^u  is  the  scalar 

T 

control  torque  on  the  i-th  mode.  The  parameters  d^ ,  and  b^  are  tabu¬ 
lated  in  Table  1. 

Letting  xt  =*  [zi  z^7  ,  and  x  =  [x^,  x£,  ...  ,  j^]T  ,  (5.1)  can  be 
written  conpactly  as 

x  *  Ax  +  Bu  ,  (5.2) 


where  A  »  diag  (Aj,  Ag,  ...  ,  Ag}  with 


(5.3) 


(5.4) 


We  eh<x>so  the  output  matrix  to  lx-  of  the  Term  0  =  ding  I  r '  c ' ,  ...  , 

I  Cj 

q^}  ,  with 


Ci  =  <Cil  * 


(5.5) 


to  reflect  the  fact  that  we  penalize  only  z^'s  in  the  performance  index. 
Equation  (5.2),  together  with  the  output  equation 


y  =  Cx  , 


(5.6) 


describes  the  system  S  in  (3.1), 

We  new  apply  to  S  the  model  reduction  scheme  based  on  the  balanced 
representation  of  S  .  However,  before  doing  this  we  note  from  Table  1  that 

and  in 

(5.3)  are  not  stable.  This  prevents  the  model  reduction  scheme  from  being 

applied  directly,  for  it  requires  the  system  be  stable.  To  avoid  this 

T  T  T 

difficulty,  we  partition  the  state  vector  as  x  =  [x^  x^]  ,  where  xy  * 

[x^  Xg  Xgl  and  xd  *  [x^  Xg  ...  x^]  correspond  respectively  to  the  un- 

danped  and  danped  modes.  Partitioning  A  and  B  of  (5,2)  and  C  of 
(5.6)  accordingly,  we  represent  S  as  the  parallel  connection  of  two  sys¬ 
tem;  S  and  S  .  ,  which  are  described  bv 
u  d 


the  first  three  modes  are  undanped,  i.e. ,  the  blocks  A^ ,  Ag 


(5.7) 


Adxd  +  BdU 


yd  =  Cd*d 


(5.8) 


Now  since  is  stable,  we  can  obtain  a  reduced  order  model  3d  for 

based  on  its  balanced  representation.  The  parallel  connection  of  and 

3d  then  gives  a  reduced  order  model  3  of  S  .  In  ether  words ,  we  keep 
the  undanped  modes  in  the  reduced  order  model  and  add  to  these  a  part  of  the 
danped  modes. 

We  choose  R  =  10~^I  ,  and  c.^  in  (5.5)  as 


c 


il 


1  i  =  1,  2,  3 

0.1  ,  i  =  4,  5,  ...  ,  9 
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to  penalize  the  undamped  modes  more  than  the  dnnped  ones.  The  optimal  gain 
matrix  and  the  eigenvalues  of  the  <pi  imal  system  are  shewn  in  Table  2,  along 
with  the  singular  values  of  the  product  EqZ ^  of  the  observability  and  con¬ 
trollability  granmians  of  Sd  .  From  these  values  we  observe  that  the  last 
four  singular  values  are  very  small  conpared  to  the  first  eight,  which  sug¬ 
gests  that  an  8-th  order  reduced  model  should  adequately  represent  Sd  .  The 
results  of  subcptimal  design  based  on  an  8-th  order  reduced  model  of  Sd  are 
shown  in  Table  3.  The  subcptimal  gains  and  eigenvalues  of  the  subcptimal  sys¬ 
tem  coincide  with  those  of  the  optimal  ones  up  to  the  fourth  decimal  place, 
producing  a  degree  of  siboptimality  u  *  1.000  and  thus  confirming  our  obser¬ 
vation.  The  results  of  the  subcptimal  design  based  on  a  second  order  reduced 
model  of  Sd  are  shown  in  Table  4,  which  shows  a,  67%  less  of  performance 
quality  corresponding  to  an  83%  reduction  in  the  order  of  S .  . 


4.6  CONCLUSIONS 


All  model-reduct ion  schemes  available  in  the  literature  are  concerned 
with  the  open-loop  behavior  of  systems,  and  therefore,  lack  the  justifi¬ 
cation  of  producing  lcwer  order  models  which  would  also  represent  adequately 
the  closed-loop  performance  of  the  systems  they  approximate.  The  suboptimal¬ 
ity  concept  outlined  in  this  paper  provides  a  means  for  validation  of  reduced- 
order  models  by  ccnparing  the  performance  of  closed- loop  reduced  model  with 
that  of  the  full-order  system.  The  advantages  of  this  validation  approach 
are: 

(i)  It  is  particularly  .suitable  Tor  measuring  effectiveness  of  con¬ 
trols  under  information  structure  constraints; 

(ii)  It  can  be  successfully  applied  to  interconnected  systems  to  pro¬ 
duce  subopt imal  control; 

(iii)  It  provides  robustness  measures  such  as  gain  and  phase  margin, 
explicit ly  in  terms  of  the  degree  of  suboptimality. 

A  drawback  of  the  suboptimality  approach  is  the  necessity  to  solve  the 
control  problen  for  the  full-order  system  to  be  used  as  reference.  This  may 
not  be  very  important  unless  the  order  of  the  systan  is  extremely  large  or 
it  contains  uncertian  parameters,  but  becomes  a  challenging  difficulty  in 
these  cases. 
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Table  1.  Parameters  of  the  model  of  the  large  space  structure 


i 

d.C  x  10"2) 

(ti. 

l 

bT(  x  10"4) 

1 

0. 

0. 

3.460 

3.460 

2 

0. 

0.148 

-9.120 

34.400 

3 

0. 

0.288 

-1.300 

0.548 

4 

1.414 

1.414 

-1.400 

36.000 

5 

1.829 

1.829 

-4.070 

-34.700 

6 

2.846 

2.846 

28.500 

39.500 

7 

3.831 

3.831 

71.800 

-24.600 

8 

8.798 

8.798 

-8.000 

x  10"4  -1.150  x 

10"5 

9 

. 

9.613 

9.613 

-0.172 

-4.130  x 

IQ'3 

Table  3.  Subcptimal  Design  I 
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Table  4.  Suboptimal  resign  II 
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5.0  PARAMETER  PLANE 


The  utilization  of  the  parameter  olane  theory  provides 

the  control  system  designer  with  an  additional  degree  of  freedom 

to  be  utilized  in  designing  control  systems.  In  a  sense  it 

is  a  two  demensional  root-locus  in  that  one  may  design  the 

control  system  through  pole  placement  techniques.  The  Control 

Dynamics  Company's  standard  practice  is  to  design  control  systems 

utilizing  the  power  of  the  parameter  plane.  Proper  use  of  the 

information  will  provide  the  control  system  desioner  with  not 
only  expected  performance  at  the  operating  point  but  one  can  also 

obtain  a  feel  of  the  relative  "robustness”  of  the  design  if 
the  selected  point  is  comfortably  away  from  the  stability 

boundari es . 

We  will  present  first  (Section  5.1)  an  overview  of  the 
theory  of  the  parameter  plane  from  the  standpoint  of  a  continuous 
system.  Then,  in  Section  5.2  the  extension  of  this  fundamental 
theory  to  include  systems  of  higher  order  will  be  covered.  Next, 
(Section  5.3)  we  will  disucss  an  overview  of  the  extension  of 
the  continuous  domain  theory  to  include  digital  systems.  Finally 
in  Section  5.4  we  will  present  some  typical  results  from  a  digital 
appl i cation . 
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5.1  OVERVIEW 

We  begin  our  study  of  the  theory  of  the  parameter  plane  by 
investigating  the  specification  of  the  stability  boundaries. 

We  see  that  there  are  three  boundaries  present,  just  as  there 
are  in  the  root  locus  theory.  There  is  a  boundary  at  s  =  0, 
a  boundary  at  ®  and  the  complex  conjugate  boundary  at  s  *  iw 
(the  imaginary  axis  in  the  root  locus).  One  represents  the 
characteristic  equation  in  terms  of  two  arbitrary  constants, 
each  of  which  multiplies  a  polonominal  in  s,  and  the  remainder 
of  the  characteristic  equation.  These  stability  boundaries  can 
then  be  represented  in  terms  of  these  two  constants  and  plotted 
as  a  function  of  these  two  constants,  hence  the  parameter  plane. 

This  is  summarized  on  Fig.  5.1. 

Now  that  one  has  the  solution  of  the  two  equations  arising 
from  setting  the  real  and  imaginary  parts  of  the  characteristic 
equation  to  zero,  one  is  made  aware  of  another  (fourth)  stability 
boundary.  This  boundary  occurs  if  the  Jacobian  is  equal  to  zero. 

In  addition  the  sign  of  the  Jacobian  provides  one  information 
on  how  the  roots  of  the  characteri s ti c  equation  behave  as  the  part¬ 
icular  boundary  is  crossed.  (See  Fig.  5.2.)  Once  the  stability 
boundaries  are  established,  then  one  is  concerned  about  the  per¬ 
formance  of  the  system.  This  performance  can  be  determined  or 
specified  by  plotting  the  contours  of  closed  loop  damping  and 
frequency  on  the  same  parameter  plane.  One  first  must  substitute 
for  the  real  and  complex  oarts  of  the  variables,  as  shown  in  Fig.  5.3. 


.  l.D-  DECOMPOSITION 


R*  j  A  (»o> )  |  -  A ,  )  «  ♦  6 ,  {u2)r  *C ,  («*)  •  0 

lm  J  A  (**>)( *Aj (ca^)  ■ 

SOLVE  FOR  «,yi 


PARAMETER  PLANE 


'=X„+iY„*l  +  iO 


Now  we  must  go  through  exactly  the  same  procedure  that  we 


did  to  determine  the  stability  boundaries.  That  is,  we  must 
equate  the  real  and  imaginary  parts  to  zero.  Thereby  we  obtain 
two  equations  in  two  unknowns  that  can  be  solved.  The  solution 
of  these  equations  then  has  a  similar  form  to  the  one  that  we 
experienced  before.  (See  Fig.  5.4.). 

Of  course,  one  must  allow  for  the  occurence  of  real  roots 
in  the  stable  region.  One  assumes  a  form  of  the  real  root  and 
makes  substitution  into  the  characteri sti c  equation.  The  result 
is  an  equation  In  two  unknowns  that  can  be  solved  parametri cal ly 
and  then  plotted  on  the  same  parameter  plane.  (See  Fig.  5.5.). 

Thus  we  have  all  the  tools  necessary  to  completely  describe 
the  roots  of  the  governing  characteri sti c  equation,  and  thus  the 
system  response,  in  terms  of  two  unknown  parameters.  We  now  can 
select  these  two  parameters  for  the  particular  system  to  achieve 
the  desired  stability  and  response  characteristics.  We  will 
wait  until  we  have  completed  the  discussion  of  the  digital  para¬ 
meter  plane  before  presenting  any  numerical  examples  of  the 
technique . 
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5.2  MATRIX  FORMULATION  EXTENSION 


5.2.1  INTRODUCTION 


In  this  section  we  focus  our  attention  on  the  mathematical 


model 

Mr|  +  Cfl+Knaf  (1) 
where  M,C,  and  K  are  nxn  matrices,  and  n  and  f  are  n-dimensional 
vectors.  The  feedback  f  is  assumed  to  be  of  the  form 


f  =  $F0 


(2) 


where  «  3  {<t>l ,  <t>z} 

"PID"  form 


0  3 


is  composed  of  two  n-vectors,  and  has  the 


(V  n) 

/  ( 4> j .  n)dt 
//  ( j .  rjdtj  dt2 

U2*  M 

/  ( 4*2  ’  n )  d  t 

//  ( 4> 2  *  n)dtj  dt2 


(3) 


and  F  =  (f.-.-)  is  a  2x6  feedback  matrix  of  constant  gains.  The 

J 

notation  (,  )  stands  for  the  inner  product  of  two  vectors. 

The  purpose  of  this  portion  of  the  report  is  to  find  a 
convenient  method  for  performing  a  parameter  plane  analysis  when 
the  parameters  of  interest  are  any  two  elements  of  the  F  matrix. 
Thus,  we  implicitly  assume  that  all  of  the  pertinent  elements  in 
the  above  three  equations  are  numerically  know  from  a  separate 
design  analysis.  One  benefit  of  the  parameter  Diane  technique 
would  then  be  a  measure  of  robustness  for  the  "design  point". 


The  analysis  that  follows  does  not  assume  any  special  form 


of  the  matrices  other  than  that  previously  described.  Thus, 
this  technique  is  itself  robust  as  it  can  be  applied  to  any  model 
which  can  be  described  by  the  equations.  The  FORTRAN  program 
developed  for  this  purpose  is  included  under  the  name  DARPA.  It 
is  very  interactive  and  only  requires  the  user  to  specify  the 
two  elements  in  the  F  matrix  and  the  range  of  values  of  the  com¬ 
plex  number  which  is  a  third  parameter  in  the  parameter  plane 
technique.  Since  the  particular  form  of  the  mathematical  model 
is  almost  diagonal,  a  simple  check  case  can  be  performed  by 
assuming  the  F  matrix  to  be  identically  zero.  Of  course,  this 
would  not  be  the  case  in  practice,  but  was  used  to  verify  the 
results  of  the  more  general  theory.  This  simplified  version  Is 
Included  as  the  FORTRAN  program  LDARPA  (Little  DARPA). 

A  second  use  of  the  parameter  plane  technique  is  to  gen¬ 
erate  two  elements  in  the  feedback  matrix  F  that  provide  a 
"design  point".  If  contours  are  plotted  for  a  fixed  damping 
factor  and  if  one  selects  a  point  of  crossing  in  the  stable  region, 
then  this  design  point  will  furnish  a  system  that  is  stable  and 
at  least  two  modes  will  have  a  known  damping  factor  and  frequency. 

An  attempt  is  then  made  to  focus  attention  on  two  parameters 
in  the  vector.  These  generalized  locations  and  slopes  are  thought 
to  be  critical  and  variations  about  the  assumed  values  are  deemed 
important.  Initial  efforts  are  centered  on  the  first  four  com¬ 
ponents  in  these  vectors,  and  (as  later  discovered)  due  to  the 
degenerate  form  of  the  sample  mathematical  model,  no  positive  results 


'I 


can  be  obtained.  Lack  of  time  prevented  further  analysis  on  this 
portion  of  the  research  but  it  is  felt  that  this  tyDe  of  analysis 
can  still  prove  to  be  fruitful  if  damping  is  included  in  the  third 
and/or  fourth  modes  or  if  components  greater  then  four  are  invest¬ 
igated. 


5.2.2  MATHEMATICAL  ANALYSIS 

The  mathematical  analysis  that  follows  assumes  the  two 
elements  of  interest  are  elements  of  the  F  matrix.  These  two 
elements  will  be  denoted  as  (at  .  8)  will  be  allowed  to  vary,  and 
the  remaining  elements  in  Eqs .  1-3  will  be  assumed  fixed.  The 

central  concept  of  the  parameter  plane  technique  is  to  discover 
those  values  of  (a, 6)  which,  for  a  fixed  complex  number,  satisfy 
the  characteristic  equation  of  the  system.  However,  the  deriv¬ 
ation  of  the  characteristic  equation  is  a  formidable  task  and  the 
alternate  technique  of  evaluating  the  determinant  of  the  Laplace 
transform  of  the  system  equations  will  yield  identical  results. 
Since  the  parameters  (a, 8)  are  spread  throughout  the  equations, 
one  must  first  determine  the  generic  (albeit  theoretical)  form 
of  the  characteristic  equation.  For  this  purpose  we  decompose 
F  as  follows:  _ 


f  f  *  f 

~ 1 2  "  1 3  1  14 

f  15  ' 

f16 

f21 

f  22  f  23  1  f  24 

f  25 

f  26 

F 


(4a) 


Bank  1 


Bank  2 


(4b) 


Bank  3 


Bank  4 


Thus,  to  say  that  a  is  in  Bank  4  will  mean  that  a  *  f  ,  or  f^g, 
or  f26. 

DEFINITION!  The  elements  a  and  $  are  in  diagonal  banks  if; 

a  is  an  element  of  Bank  1  and  6  is  an  element  of  Bank  4, 
or  8  is  an  element  of  Bank  1  and  a  is  an  element  of  Bank  4, 

or  a  is  an  element  of  Bank  2  and  8  is  an  element  of  Bank  3, 

or  8  is  an  element  of  Bank  2  and  a  is  an  element  of  Bank  3. 

Otherwise,  a  and  6  are  in  nondiagonal  banks. 

The  generic  form  of  the  characteristic  equation  assumes  one 
of  two  forms  and  the  choice  of  the  form  is  dependent  on  the  location 
of  a  and  6.  More  specifically,  we  have  the  following: 

THEORM  (i)  If  (  a, 8)  are  in  nondiagonal  banks,  then  the  charact¬ 
eristic  equation  of  Eqs.  1-3  is  equivalent  to 

A(s)  *  1  +  ab(s)  +  Bc(s)  *  0,  (5) 

where  b(s)  and  c(s)  are  rational  polynomials  Independent  of 


(a, 6) . 


(ii)  If  (a, 8)  are  in  diagonal  banks,  then  the  characteristic 


equation  is  equivalent  to 

A ( s )  =  1  +  ab(s)  +  8c(s)  +  a6d(s)  (6) 

where  b(s),  c(s),  and  d(s)  are  rational  polynomials  independent  of 

(a,8) . 

REMARK.  The  oarameter  plane  loci  is  now  determined  by  selecting 
convenient  values  of  (a, 8)  to  evaluate  the  left  side  of  Ea.  5  or 


Eq.  6.  The  information  is  then  used  to  find  b(s),  c(s),  and/or 
d(s),  from  which  (a,S)  is  determined  by  solving  Eq.  5  or  Eq.  6. 

This  portion  of  the  analysis  is  continued  in  the  next  section. 

The  reader  not  interested  in  the  detailed  mathematical  analysis 
of  proving  the  above  theorem  may  wish  to  skip  to  that  section. 

The  proof  of  the  theorem  depends  on  the  following  well 
known  lemma: 

LEMMA  1.  Let  I  denote  the  identity  matrix  and  b,c,  be  nonzero 
vectors.  Then 

det( I  -  bcT)  =  1  -  (c,b)  (7) 

PROOF.  Since  (I  -  bcT)b  =  {1  -  (c,b)}b,  we  have  that  b  is  an 
eigenvector  of  I  -  bcT  with  eigenvalue  1  -  (c,b).  Let  x  be  any 
other  eigenvector  with  corresponding  eigenvalue  s,  so  that  (I- 
bcT)x  *  sx  implies  (1  -  s)x  »  (c,x)b.  Thus,  if  s  is  not  equal 
to  1  then  x  may  be  selected  as  b  above.  If  s  -  1,  then  the 
equation  (c,x)  *  0  has  n-1  independent  solutions.  Thus,  the 
spectrum  of  I  -  bcT  consists  of  the  two  eigenvalues  1  -  (c,b) 
and  1.  The  proof  of  the  lemma  is  completed  by  noting  that  the 
product  of  the  eigenvalues  of  a  matrix  is  equal  to  the  determinant 
of  the  matrix. 

We  also  require  the  following  extension  of  Lemma  1. 

LEMMA  2.  Let 

C  =  I  -  bcT  -  deT 

and  assume  the  vectors  b  and  d  to  be  independent  as  well  as  c  and 
e.  Further  assume  that  (c,d)  t  0  and  (e,b)  f  0.'  Then 

det(C)  =  I  -  (c.b)  -  ( e ,  d )  +(c,b)(e,d)  -  (c,d)(e,d).  (8) 
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PROOF.  As  in  the  proof  of  Lemma  1,  we  show  that  the  spectrum  of 
C  is  composed  of  the  three  eigenvalues  {l.s^.sj  where 

s .  *  1  -  (c,b)  -  u  .  (c,d) ,  j  *  1.2 

J  J 

and  ui  are  roots  of  the  quadratic  equation. 

(c,d)u2  +  { (c ,d)  -  (e,d)  }u  -  (e,b)  «  0.  C9) 

For  this  purpose,  we  observe  that  there  are  n-2  independent 
eigenvectors  with  eigenvalues  eaual  to  one.  That  Is,  Cx  =  x 
yields 

x  -  b(c,x)  -  d(e,x)  =  x  or 
b(c,x)  +  d(e,x)  *  0. 

Since  b  and  d  are  independnet  vectors,  we  conclude  from  this 
last  equation  that  (c,x)  =  0  and  (e,x)  «  0.  The  assumed  inde¬ 
pendence  of  c  and  e  Implies  that  there  are  n-2  distinct  solutions 
to  these  two  last  equations. 

Now  let  x  =  b  +  ud  so  that  Cx  *  sx  and  the  independence 
of  b  and  d  yield 

1  -  (c,d).-  u(c,d)  =  s,  and 
u  -  (e,b)  -  u(e,d)  =  su. 

If  we  multiply  the  first  equation  by  -u  and  add  to  the  second 
equation,  we  obtain  Eq.  9.  Now  the  roots  of  Eq.  9  satisfy 
-(c.dMuj  +  u?)  =  (c,b)  -  (e,d),  and 

(c.dJUjUg  =  - (e ,b ) . 

Thus , 

det(C)  *  s j s ^  *  [(  1  -  (c,b)}  -  u^c.djJQl  -  (c,b)} 

-  u2  (c.dfl. 
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and  Eq.  8  follows  by  expanding  this  last  result. 

REMARK.  The  assumption  of  the  vectors  being  nonorthogonal  can  be 
relaxed,  but  this  is  not  required  in  this  report. 

We  now  return  to  the  proof  of  the  Theorem. 

PROOF,  (i)  To  solidify  ideas  let  a  3  f  and  6  =  f^  so  that 
(a, 6)  are  in  Bank  1.  Then 

r~  H 

a  6  0  0  0  0 

F  3  F'  + 

0  0  0  0  0  0 

where  the  2x6  F'  matrix  is  formed  from  F  bv  deletions  a  and  6  . 

Thus, 

Fe  =  F'G  +  {0(44,  n)  +  S/U,  ,n)dt  0}T 

and 

f  3  $F'0  +  ad>1  (d>!  ,n)  +  8*T/Ultn)dt 
To  simplify  notation  we  define  the  nxn  matrix  R(s)  as  follows: 

R(s)  3  Ms2  +  Cs  +  K  -  *F'  {sd>  v s-1^*  <t>2,  <t>z,  S-1^}. 

Then  the  Laplace  transform  of  Eqs.  1-3  yields. 

{ R  ( s  )  -  (as  +  3)  4»[  4>  1 T }  n  (s)  3  I.V.  (10) 

where  n(s)  is  the  Laplace  transform  of  n(t)  and  I.V.  stands  for 
initial  values.  Thus,  the  characteristic  polynomial  of  the  system 
becomes 

p(s)  3  det  rR(s)  -  (as  +  8)  4^^}. 

If  we  assume  the  inverse  of  R(s)  exists,  then  we  can  factor  out 
R(s)  on  the  left  to  obtain  the  identity 

R(s )  -  (as  +  S)*!*!1  =  R( s }  (I  -  (as  +  3} R'  1  ( s ) 4 , 4  l  T} . 
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Usinn  Lemma  1  with  b  *  (as  +  6)  R* 1  ( s  )<j>  l  ,  c  =  <J>  ,  and  using 

the  last  identity  yields  the  following  form  for  the  characteri s ti c 

polynomial 

P  ( s  )  *  det  { R  ( s  ) }  {1  -  (as  +  8)  (d>1  ,R-1(s)4>i)>. 

Thus,  p ( s )  »  9  implies  the  characteri sti c  equation  can  be  expressed 
as 

A  ( s )  *  1  -  (s  (4>  t » R  1  ( s )  ! )  >o  -  (d>1,R“1(s)d>I)6. 

In  general,  the  inverse  of  R { s )  exists  as  it  is  composed  of  rational 
polynomials,  and  the  theorem  is  proved  for  the  special,  case  where 
c(s)  *  -($i.R  1  ( s )  i )  and  b(s)  B  sc(s).  The  remaining  parts  of 
(i)  can  be  proved  by  moving  the  elements  a  and  8  to  their  specified 
locations  and  repeating  this  proof. 

(1i)  Again  to  solidify  ideas,  we  let  a  =  f^  be  in  Bank  1 
and  6  *  be  in  Bank  4.  Then 

( <t> ,  »n)a 

f  =  *F'0  +  * 

(<fr2  n)B 

where  F'  is  formed  from  F  by  deletion  a  and  8.  Thus 

P ( s  )  =  det  { R ( s  )  ’  -  as<j>  1  $  j T  -  6s<j>2<P2T> . 

Using  Lemma  2  with  b  =  -  asR'1^,  c  =  <p t ,  d  =  -SsR"1^ 

and  e  *  <t>2  yields 

p ( s )  *  det  { R ( s )  HI  -  as(^j  ,R'1(t>i )  -  Ss(d>2,R-1<J>2) 

+  a8  G($  j » R* ^ t )  (*t>2  , R"  1d>2 )  -  (<f>i  , R”  1d>2 )  (d>2  ,R_1<1),  3>s2 

so  that 

b(s  )  *  -s  ( d>  x  ,R_10i), 
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c ( s )  =  -s(4>2,R_14>2},  and 

d(s)  3  s2{(<f>1  ,R"1<|>1)(d>2,R'1<i>2)  -  (*i.R"1«2)(*a,R’1^x)}. 

5.2.3  COMPUTER  PROGRAMS 

In  general,  the  R  matrix  as  defined  by  Eq.  10  will  not 
be  sparse  and  the  evaluation  of  b(s),  c(s),  and  d(s)  in  the 
Theorem  will  require  considerable  computational  effort  to  determine 
its  inverse.  To  alleviate  computing  this  inverse  for  each  value 
of  s,  we  return  to  the  original  form  of  the  characteristic  poly¬ 
nomial.  That  is, 

p(s  ;o,B)  3  det{Ms2  +  Cs  +  K  -  $F(ct ,  8)0( s ) } .  (11) 

Then  using  Eq.  11  with  Eq.  5  (or  Eq.  6)  yields 
p(s;0,0)  3  det  {R(s) } 
p(s  ;  1 ,0)  *  det  { R ( s  )  Hi  +  b ( s  ) } 

p  ( s  ;  0 , 1 )  *  det  {R(s )  Hi  +  c(s)}  (12) 

(and  for  Eq.  6:  p(s;l,l)  3  det  '(R(s)Kl  +  b(s)  +  c(s)  +  d(s)}). 

The  information  on  the  left  of  Eq.  12  is  found  by  numerically 
evaluating  the  complex  determinants  by  using  a  modified  Gaussian 
elimination  scheme  to  handle  complex  elements.  As  Eq.  5  or  Eq.  6 
gives  the  explicit  form  of  the  characteri s ti c  equation,  we  can  now 
solve  for  b(s),  c(s),  and,  if  needed,  d(s).  Finally,  the  values 
of  a  and  8  that  satisfy  the  characteri sti c  equation  are  found 
by  taking  the  real  and  imaginery  parts  of  Eq.  5  or  Eq.  6  and 
solving  the  resulting  equations.  For  Eq.  5,  the  resulting  equations 
are  linear  in  the  two  unknowns  while  for  Eq.  6  the  equations  are 
nonlinear  and  require  a  more  delicate  routine  for  their  solutions. 

These  subroutines  have  I een  denoted  as  S0LVE1  and  S0LVE2 ,  respectively. 
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The  analysis  of  the  problem  has  been  coded  in  FORTRAN 
using  double  precision  arithmetic.  Extensive  use  is  made  of 
complex  arithmetic,  but  for  the  sake  of  program  portability  no 
use  is  made  of  FORTRAN  complex  variables.  The  main  line  program 
is  name  DARPA  and  follows  the  presentation  of  the  previous  para¬ 
graph.  DARPA  requires  the  initial  matrices  M,C,K,<&,  and  F 
from  the  user's  data  base  and  then  requests  for  the  row  and 
column  of  the  desired  elements  in  the  F  matrix  that  are  to  be 
used  in  the  parameter  plane  analysis.  Then  the  user  is  prompted 
for  a  range  of  values  of  the  variable  s.  DARPA  calls  subroutines 
BUILD,  MMULTD ,  CDIV,  COND,  CGAUSS,  S0LVE1,  AND  S0LVE2 .  Subroutine 
BUILD  is  used  to  construct  the  real  and  imaginary  parts  of  the 
matrix  in  Eq.  11  with  F  replaced  by  F'.  (The  present  version  of 
BUILD  assumes  the  special  form  of  the  mathematical  model  for 
test  purposes.)  MMUlTD  is  a  double  precision  multiplication 
routine,  DCIV  is  a  double  precision  comolex  variable  division 
routine,  CGAUSS  is  a  modification  of  Gaussian  elimination  to 
handle  complex  numbers.  For  some  values  of  s,  round  off-error  may 
be  experienced,  so  a  conditioning  routine  COND  is  included  to 
scale  the  elements  of  the  matrices  before  calling  GAUSS.  Finally, 
S0LVE1  and  S0LVE2  as  mentioned  in  the  previous  paragraph  are  used 
to  solve  the  resulting  equations. 

A  sample  run  of  DARPA  is  included  after  the  listings  of  the 
program.  The  first  part  of  the  run  lists  the  current  numerical 
values  of  the  input  matrices  (F  was  assumed  to  be  zero  for  the 
test  case).  The  second  part  of  the  output  prompts  the  user  for 
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the  values  of  a  and  S.  In  the  example  a  *  and  S  *  f  were 
selected  as  the  parameters  of  interest.  Then  the  user  is  prompted 
for  initial  values  of  s  =  SIG  +  j  OMEGA,  incremental  values  of 
s  and  final  values  of  s.  In  the  example  at  hand  two  points 
s  -  .lj  and  s  -  .2j  were  run.  Note  in  the  example  that  the 
increment  in  SIG  was  zero  so  a  nonzero  value  was  selected  for 
the  final  value  of  SIG.  The  results  of  the  calculations  are 
given  next.  In  order  to  determine  some  accuracy  of  results,  the 
calculated  values  were  substituted  back  into  the  original  equation 
and  the  determinant  was  evaluated  once  aaain.  This  determinant 
should  be  zero.  The  final  calculation  shows  this  diagonal  elements 
of  this  determinant  after  calling  CGAUSS.  As  the  determinant  will 
be  the  product  of  these  diagonal  elements,  it  is  seen  that  the 
results  are  within  predictable  accuracy  (note  the  last  two  elements) 
To  gain  further  confidence  in  the  DARPA  code  a  separate  and  inde¬ 
pendent  analysis  was  performed  on  the  specific  form  of  the  equations 
This  is  described  in  the  next  few  paragraphs,  and  by  comparing  the 
results  of  DARPA  and  LDARPA  we  note  the  excellent  agreement  of  the 
numerical  results. 

If  we  assume  F'  -  0  then  the  R(s)  matrix  above  Eq.  10  is 
R(s)  =  Ms2  +  Cs  +  K. 

As  R( s )  for  the  test  case  is  almost  diagonal  an  explicit  expression 
for  its  inverse  is  easy  to  obtain  and  one  can  obtain  a  much  simplier 
form  of  the  characteristic  equation  of  the  closed  loop  system.  For 
example,  if  a  =  f  and  6  =  f  we  have  already  seen  that 
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A(s)  *  1  -a{sU,,R  ^s)#,)}  -6(4>1,R'*1(s)<t»1). 

Substituting  R“ 1 ( s )  in  this  last  expression  and  expanding 
the  inner  products  gives  us  the  following  expression  for  the 
characteristic  equation: 

10  2  r  2  1 

A ( s )  =  1  -  (as  +  6)  {?_  <*>!  (j)  /  [s  +  c(j)s  +  k(j)J  (13) 

J  *  ^ 

+  Cmd.l)  +  m(2,2)  -2m(l,2)]  /  [s2  (m(l,l)m(2,2)  -  m(l,2)2Q  } 

where  c(j)  and  k(j)  are  the  diagonal  elements  of  C  and  K,  respectively. 
This  formulation  allows  for  a  much  smaller  program  with  much  less 
chance  for  numerical  roundoff.  As  mentioned  above  the  coding  of 
Eq.  13  is  called  LDARPA . 

5.2.4  DISCUSSION 

Fig.  1  is  a  sketch  of  the  parameter  plane  with  a  =  o.  The 

parameters  selected  were  a  *  f,,  and  B  =  f  .  The  three  Dages  with 

11  12 

the  title  PRODUCTION  RUN  FOR  SIGMA  *  0  are  the  actual  numerical 
values  used  in  producing  this  sketch  and  were  generated  using 
LDARPA.  Note  that  this  sketch  is  not  to  scale,  due  to  the  large 
dispersion  of  the  magnitude  of  the  numerical  values.  This  sketch 
separates  the  plane  into  stable  and  unstable  regions.  The  static 
boundary  (the  f  -  axis)  was  determined  by  setting  s  *  0  and 
performing  a  separate  analysis.  The  stable  region  was  determined 
by  reducing  the  problem  to  one  of  smaller  dimension  and  invoking 
the  Routh-Hurwi tz  stability  test.  For  all  intents  the  stable 
region  consists  of  the  third  quadrant.  Fig.  2  is  an  enlargement 
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of  the  region  close  to  the  origin.  It  merely  shows  that  for  small 

positive  values  of  f  the  system  will  remain  stable  provided 

f  is  negative. 

12 

Also  included  under  the  title  PARTIAL  PRODUCTION  RUN  FOR 
SIGMA  *  -.01414  are  some  numerical  values  for  the  parameter  plane 
with  a  fixed  and  equal  to  the  value  in  the  title.  Lack  of  time 
and  space  prohibits  exhibiting  the  graph  of  this  data  and  the 
other  numerious  numerical  values  generated.  However,  it  will  be 
noted  that  a  careful  scan  was  made  of  the  data  for  a  *  -.707  and 
no  multiple  crossing  of  the  graphs  were  detected. 

Some  attempt  was  made  to  perform  a  parameter  plane  on  the 
elements  in  the  $  matrix.  For  a  full  F  matrix  the  resulting 
equations  appear  to  be  highly  nonlinear  if  two  arbitray  values 
are  selected,  hence  no  numerical  experiments  were  performed.  How- 

t 

ever,  for  the  special  model  considered,  and  if  only  two  values  of 

the  F  matrix  are  nonzero  then  one  can  still  get  some  partial  results. 

For  this  purpose,  let  us  assume  f  and  f  are  the  nonzero  elements 

11  12 

and  and  concentrate  on  the  first  bending  mode.  That  is  let 


x  *  4>j(3)  and  y  *  <t>2(4) , 

so  that  Eq.  13  can  be  rewritten  as 


s2  +  c ( 3 ) s  +  k( 3) 

m(  1  a),  +_m(2_,2)  _-_2mj  1 ,2) 
m (l,l)m(2,2)  -  m(l,2)2 


s2  =  c (4 ) s  +  k ( 4 )  (f  s  +  f  ) 

11  12 

.  ;• . ?l(ii . 

j=5  s2  +  c(j)s  +  k(j)  •. 
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From  the  given  numerical  data  c(3)  *  c(4)  *  0.  If  one  now  sub¬ 
stitutes  s  *  jw  in  this  last  expression  we  run  into  a  contradiction 
as  the' imaginary  part  of  the  right  hand  side  must  be  zero.  To 
alleviate  this  problem  one  must  assign  appropriate  values  to  the 
modal  damping  matrix  C. 


79 


o  o  o  o  o  -L 


DARPA  LISTING 


DOUBLE  PRECISION  MASS < 4) . C< 6) , STIFF ( 3 ) ,PHI1 ( 10) ,PHI2< 10) 

DOUBLE  PRECISION  F(2, 6) . X ( 10, 10) , Y< 10, 10) , WORM 10, 10) ,U<6, 10) 
DOUBLE  PRECISION  SIG, W, DSIG, DW, FSIG, FW, AR, AI , BR, BI . CR. Cl , DR, DI 
DOUBLE  PRECISION  ALPHA, BETA, ALPHA 1 , ALPHA2, BETA 1 , BETA2 
DOUBLE  PRECISION  MAXA, MAXB, MAXC, MAXD, RATIO 
DIMENSION  I PVT (10) 

READ (6, 100)C 
READ <6, 100) ST IFF 
READ (6, 100) PHI 1 
READ <6, 100) PH I 2 
READ (6, 100) MASS 
N  =  10 

DO  4  I  =  1,N 
DO  4  J  =  1 ,  N 
F  <  I , J )  =  O.ODO 
CONTINUE 

...CURRENTLY  ASSUMING  F  IS  ZERO.  IF  NOT,  UPDATE  DATABASE  AND  INSERT  A 
. . .READ <6, 100)  F.  NOTE  THAT  MASS, C, AND  STIFF  ARE  VERY  SPECIAL  AND 
...ONLY  INCLUDE  THE  NONTRIVIAL  ELEMENTS.  IF  THEY  WERE  GENERAL  THEN 
...WILL  NEED  TO  CHANGE  THE  INPUT  OUTPUT  PORTIONS  OF  THE  PROGRAM. 
FOLLOWING  OUTPUT  DATA  FORMAT  IS  FOR  THIS  SPECIAL  CASE 

100  FORMAT ( 10D11.6) 

WRITE (5, 200) MASS 

200  FORMAT  <  "  Ml  1 ,  M12 ,  M21 ,  M22  *  ',/,4D12.6) 

WRITE (5, 201 >C 

201  FORMAT ( "  LAST  6  ELEMENTS  OF  DIAG.  DAMPING  MATRIX  =  ',/,6D12.6> 
WRITE (5, 202) 

202  FORMAT ('  LAST  S  ELEMENS  OF  DIAG.  STIFF.  MATRIX  *  ') 

WRITE <5, 203) ST IFF 

203  FORMAT <5D 12. 6) 

WRITE  <5, 204 ) 

204  FORMAT ( "  PH II  =  '  ) 

WRITE (5, 203) PHI 1 
WRITE <5, 205) 

205  FORMAT  (  •'  PH  1 2  =  '  ) 

WR I TE ( 5 , 203 ) PH I 2 
WR I TE ( 5 , 206 ) 

206  FORMAT (//'  CURRENT  RUN  HAS  F  =  0  ') 

WRITE (5, 207 ) 

207  FORMAT ( / '  INPUT  THE  ROW  AND  COLUMN  OF  ALPHA' / ) 

READ  (4,208) I, J 

208  FORMAT (2 12) 

WRITE <5, 209) 

209  FORMAT (  /'  INPUT  THE  ROW  AND  COLUMN  OF  BETA-/) 

READ  (4,203)  K, L 

C 


DETERMINE  IF  CA 

SE  I 

OR  CASE  I I 

I CASE  ■  1 

IF  (  (  I.  EG).  1 

.  AND. 

J.LE.3) . AND. (K.EQ. 

.AND. 

L .  GE . 4 ) ) 

I  CASE 

.  2 

IF ( ( I.EQ. 1 

.  AND. 

•J .  GE .  4  ) .  AND .  ( K .  EQ . 

a!  ■ 

AND 

.L.LE.3)  ) 

I  CASE 

SE  £ 

IF ( ( I.EQ. 2 

.  AND. 

J. LE. 3) . AND. <  K . EQ. 

1  . 

AND 

.  L.0E.4)  ) 

I  CASE 

-  2 

IF ( ( I.EQ. 2 

.  AND. 

•J .  uE .  4  ) .  AND .  ( K .  EQ . 

\ 

*  • 

AND 

L.LE.3)  ) 

I  CASE 

=  2 

1  WRITE (5, 300) 

300  FORMAT (/  -  ENTER  VALUES  OF  SIG  AND  OMEGA  " / ) 
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READ<  4,  301 )  SIG.W 

301  PORMAT (2D1&. S) 

WRITE <5, 302) 

READ ( 4 i 30 1 ) DS I G , DW 
WRITE <5, 303) 

302  ■  FORMAT </  '  ENTER  VALUES  OF  DELTA  SIG  AND  DELTA  OMEGA  '/) 

303  FORMAT (/  '  ENTER  VALUES  OF  FINAL  SIG  AND  FINAL  OMEGA  '/) 
READ( 4. 301 )FSIG. FW 

2  F ( I » J )  *  O.DO 

F(k'rL)  =  O.DO 

CALL  BUILD(MASS,C,STIFF,F,PHI1,PHI2,SIG,W,W0RK,U,X,Y) 

CALL  COND(N,  X,  Y.MAXA) 

CALL  CGAUSS(N,X,Y, I PVT) 

CALL  DET ( IPVT ( N) ,  N,  X  ,  Y, AR. AI ) 

F(I..J)  =  1 .  DO 

CALL  BUILD< MASS, C. STIFF, F, PHI 1.PHI2, SIG, W, WORK, U,X,Y) 

CALL  COND(N, X, Y.MAXB) 

CALL  CGAUSS <  N ,  X ,  Y ,  I PVT ) 

CALL  DET(IPVT(N),N,X,Y,BR,BI) 

F( I » J)  »  O.DO 
F(K»L)  =  1 . DO 

CALL  BUILD ( MASS, C, STIFF, F, PHI 1,PHI2, SIG, W, WORK, U,X,Y) 

CALL  COND  <  N , X , Y , MAXC  > 

CALL  CGAUSS(N,X,Y, IPVT) 

CALL  DET < I PVT <N),N,X,Y, CR , C I ) 

IF ( I CASE. EQ. 2)G0  TO  20 

CALL  CD I V  <  AR , A I , BR . B I , ALPHA 1 , ALPHA2 ) 

RATIO  *  < MAXB/MAX A ) **N 
BR  »  ALPHA 1*RAT 10  -  l.ODO 
BI  »  ALPHA2*RAT I 0 

CALL  CD  I V  ( AR ,  A I ,  C.R ,  C I ,  ALPHA  1 ,  ALPHA2  > 

RATIO  =  <  MAXC/MAXA ) **N 
CR  a  ALPHA 1*R AT 10  -  l.ODO 
Cl  =  ALPHA2*RAT I 0 
AR  a  l.ODO 
AI  =  O.ODO 

CALL  S0LVE1 <AR, AI,BR, BI,CR, Cl, ALPHA, BETA) 

C _ START  CHECK  OUT  MODE. ..CAN  BE  DELETED  IF  DESIRED. 

F ( I , J )  =  ALPHA 
F(K,L)  a  BETA 

CALL  BUILD (MASS, C, STIFF, F, PHI 1,PHI2. SIG, W, WORK, U,X,Y) 

CALL  COND <N,X,Y,MAXC) 

CALL  CGAUSS <  N , X , Y , I P VT ) 

CALL  DE T ( I P VT ( N ) , N , X , Y , AR , A I ) 

WRITE ( 4, 30 )  SIG, W, ALPHA, BETA 

30  FORMAT  ('  S I G ,  W ,  ALPHA  .BETA  =  ,  4D 1 5 .  7 ) 

WRITE (4,31) 

31  FORMAT  <  y  X(1,1),Y(1,1),X(2,2),Y(2,2)  ...  *  ' , / ) 

WR I TE ( 4 , 32 )(X(II,II),Y(II,II)»  II  =  1,N) 

32  FORMAT ( 10DS. 2) 

C...END  OF  CHECK  OUT  MODE 

WRITE  (3, 400 ) S I G , W , ALPHA , BETA 
SIG  =  SIG  +  DS I G 
W  —  W  +■  DW 

IFCSIG.LT. FSIG. AND. W.LT.FW)GO  TO  2 
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GO  TO  3 
C 

C . THIS  IS  CASE  II 

C 

20  F(I,J)  *  l.ODO 

CALL  BUILD (MASS, C. STIFF. F, PHI 1 , PHI2. SIG, W, WORK, U, X,  Y ) 

CALL  CGAUSS < N. X, Y, IPVT ) 

CALL  DET ( IPVT (N),N»X,Y,DR»DI> 

DR  *  DR  -  AR  -  BR  -  CR 
DI  »  DI  -  AI  -  BI  -  Cl 

CALL  SOLVES  <  AR , A I , BR , B I . CR , C I , DR , D I , ALPHA 1 , BETA 1 , ALPHAS , BETAS > 
WR I TE  <  8 , 500 ) S I G , W , ALPHA 1 , BET A 1 , ALPHA? , BETAS 
SIG  =  SIG  +  DSIG 
W  =  W  +  DW 

IF(SIG. LT. FSIG. AND. W.LT . FW)GO  TO  2 
3  WRITE <5, 600) 

600  FORMAT (  '  ENTER  1  TO  CONTINUE  ') 

READ ( 5, 601 )  II 

601  FORMAT (SIS) 

IF  (II.EQ. 1)  GO  TO  1 

400  FORMAT C4D1 6. S> 

500  FORMAT (6D1 6. S> 

STOP 

END 
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DARPA  SUBROUTINES 


SUBROUT I NE  BU I LD  < MASS , C , ST I FF , F , PH I 1 , PH I 2 . S I G , W , WORK , U»  X, Y) 
DOUBLE  PRECISION  MASS (4) » C<6) , STIFF<S) , PHI 1(10). PHI2( 10) 
DOUBLE  PRECISION  F<2,6) , X( 10, 10) , Y(10, 10 > , WORK( 10, 10),U<6. 10) 
DOUBLE  PRECISION  SIG, W, Cl > C2, C3. C4 
Cl  *  SIG/(SIG*SIG  +  W*W> 

C2  *  SIG*SIG  -  W#W 
C3  =*  -W/(SIG*SIG  +  W*W) 

C4  »  2. DO*SIG*W 
DO  1  I  *  1.10 
DO  1  J  =  1,10 
X<I,J)  =  O.ODO 
Y < I , J )  *  O.ODO 
CONTINUE 
DO  2  I  *  1,10 
DO  2  J  *  1.6 

WORK ( I .  J )  *  PH II  { I  >*F ( 1 ,  J>  PHI2(  I  )*F(2,  J) 

CONTINUE 

DOS  J  =  1,10 

U(1,J)  -  -SIG*PHIi C J) 

U(2, J)  *  —PHI 1 ( J) 

U(3,J)  *  — C1*PHI 1 ( J) 

U(4, J)  *  — SIG*PHI2< J) 

U<5,J)  =  -PHI2( J) 

U(6» J)  »  -C1*PHI2(J) 

CONTINUE 

CALL  MMULTD(WORK,U,X, 10,6. 10) 

X<1,1)  *  MASS <  1 ) *C2  +  X<1,1) 

XC1.2)  »  MA3S<2)*C2  .+  X(l,2) 

X (2, 1 )  =*  MASS ( 3 ) *C2  +  X(2. 1) 

X (2, 2)  »  MASS < 4 ) *C2  +  X(2.2) 

X(3,3)  «  C2  +  STIFF ( 1 )  +  X(3»3) 

X (4, 4)  »  C2  ♦  STIFF (2)  ♦  X(4,4) 

DO  5  I  »  5. 10 
KK  a  1-2 
KKK  a  1—4 

X(I,I)  a  C2  +  STIFF(KK)  +  X ( I , I )  +  SIG*C(KKK) 

CONTINUE 
DO  6  J  =  1, 10 
U(1,J)  -  — W*PH 1 1 <  J ) 

U(2, J)  *  O.ODO 
U  <  3 , J )  =  -C3*PHI1 < J) 
i  J  ( 4 , J )  =  -W*PHI2<J) 

U(5,J)  a  O.ODO 

IJ  ( 6 ,  J )  =  -C3*PHI2CJ> 

CONTINUE 

CALL  MMULTDCWORK.U, Y, 10,6,  10) 

Y ( 1 ,  1  )  =  MASS  < 1 )*C4  +  Y  < 1 , 1 ) 

Y(l, 2)  *  MASS (2)  *C4  +  Y< 1 ,2) 

Y (2, 1 )  a  MASS ( 3 ) *C4  +  Y(2,l) 

Y<  2, 2)  =  MASS  <  4 ) *C4  +  Y(2,2) 

Y  <  3, 3 )  *  C4  +  Y  <3,  3 ) 

Y ( 4, 4)  =  C4  +  Y<  4, 4 ) 

DO  I  a  5,10 

KK  =  1-4 

Y  <  I ,  I )  a  C4  +  C ( KK  >  *W  +  Y ( I , I > 
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S  CONTINUE 

RETURN 
END 


SUBROUTINE  MMULTD< A, B . C. L, M, N > 

DOUBLE  PRECISION  A<L,M> ,B(M,N> .C(L,N) 
DO  1  I  ■  1,L 
DO  1  J  a  1,N 
C(I.J)  a  O.ODO 
DO  1  K  a  1,M 

C(I»J)  =  C(I,J>  +  A< I,K)*B(K, J) 

1  CONTINUE 
RETURN 
END 


SUBROUTINE  CDIV(X,Y.U,V,AR,AI) 

DOUBLE  PRECISION  X, Y, U, V, AR, AI , D 

D  =  X*X  Y*Y 

IF  (D.EQ. O.ODO)  RETURN 

AR  *  <X*U  ♦  Y*V)/D 

AI  a  <X*V  -  Y*U)/D 

RETURN 

END 

SUBROUTINE  COND <N, X , Y, MAX ) 

DOUBLE  PRECISION  X(N,N),Y(N,N), MAX 

. IF  TROUBLE  EXISTS  WITH  CGAUSS  WE  MAY  HAVE  TO  SCALE 

....  THE  ELEMENTS  OF  THE  MATRIX  BEFORE  DOING  A  GAUSSION 
....  ELEMINATION.  OTHERWISE,  THIS  IS  A  DUMMY  ROUTINE  AND 
....  MAY  BE  IGNORED.  CURRENT  VERSION  SCALES  EACH  ELEMENT  BY  MAX. 
MAX  »  100. ODO 
DO  1  I  a  l,N 
DO  1  J  a  l.N 
X(I,J)  a  X ( I , J) /MAX 
Y( I, J)  a  Y(I,J)/MAX 
1  •  CONTINUE 

RETURN 
END 


SUBROUTINE  CGAUSS < N, X, Y, I PVT) 

DIMENSION  IPVT(N) 

DOUBLE  PRECISION  X(N,N) , Y(N,N> 

DOUBLE  PRECISION  D1 , D2, T1 , T2, T 

IF(N.EQ.l)  RETURN 

IPVT(N)  a  i 

NM1  =  N- 1 

DO  35  K  =  l,NMl 

KP1  a  K+l 

M  =  K 

DO  15  I  =  KPl,N 

01  =  X(I,K)*XCI.K)  +  Y(I,K)*Y(I.O 
02  =  X  ( PI,  K )  *X  (M,  K )  +  Y  <  M,  K)  *Y  <  M,  K) 
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IF(DABS(D1 > . GT. DABS (02) >  M  »  I 
15  CONTINUE 

IPVT(K)  *  M 

IF(M.NE.K)  I PVT (N >  *  -IPVT(N) 

T1  *  X(M,K) 

T2  =  Y(M,K) 

X(M,K)  *  X(K.K) 

Y(M,K)  *  Y(K.K) 

X (K, K)  *  T1 
Y  ( K ,  K )  *  T2 
T  =  T1*T1  +  T2*T2 
I F ( T . EQ . 0 . ODO ) GO  TO  35 
DO  20  I  *  KPl.N 

D1  *  —  (  X  ( I ,  K ) *T1  +  Y(I,K)*T2)/T 
Y ( I , K)  *  (  X  ( I ,  K ) *T2  -  Y(I,K)*T1)/T 
XCI.K)  =  D1 
20  CONTINUE 

DO  30  J  =  KPl.N 
T1  *  X(M.J) 

T2  =  Y(M. J ) 

X(M.J)  a  X (K, J) 

Y(M,J)  =  Y<K.J) 

X(K,J)  »  T1 
Y(K.J)  =  T2 
T  =•  T1*T1  +  T2*T2 
IF(T.EQ.O.ODO)  GO  TO  30 
DO  25  I  *  KPl.N 
Dl  a  X(l.K)*Tt  -  Y(I,K)*T2 
D2  *  X(I,K)*T2  +  Y< I»K)*T1 
X(I.J)  *  X<I.J>  +  Dl 
Y(I.J)  »  Y < I » J)  +  D2 
25  CONTINUE 
30  CONTINUE 
35  CONTINUE 
RETURN 
END 

SUBROUTINE  SOLVE 1 < AR. A I , BR. BI . CR, Cl , XI . X2) 

DOUBLE  PRECISION  AR, AI , BR, BI . CR, Cl . XI , X2, DET 

DET  a  0R*CI  -  BI*CR 

IF (DABS (DET) .LT. 1.0D-37)GO-  TO  10 

XI  =  -Cl /DET 

X2  »  BI/DET 

RETURN 

10  WRITE(5, 100) 

100  FORMAT ( '  TROUBLE,  LET  ALPHA  =  BETA  a  o  ) 

XI  =  0. ODO 
X2  =  0. ODO 
RETURN 
END 
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SUBROUTINE  30LVE2< AR, AI , BR, BI , OR, Cl , DR, DI , X 1 , X2> X3,  X4 ) 
DOUBLE  PRECISION  AR, AI , BR, BI , CR, Cl . DR, DI , XI , X2. X3. X4 
DOUBLE  PRECISION  A< 3) , B ( 3) , D4, BB, CC. DISC. DET, BIOA, BIOB 
A  <  3)  =  DR 
A  <  7 )  =  DI 
AO)  =  CR 
A (5 )  =  Cl 
A  ( 4 )  =*  BR 
A  <  3 )  =  BI 
A(2)  =  AR 
A( 1 )  =  AI 

A ( 3 ) XY  +  A ( 6 ) x  +  A(4) Y  +  A(2>  »  0 
A ( 7 ) XY  +  A  ( 5 )  X  +  A  <  3 )  Y  +  A<1)  =  0 

CONTINUE 

D4  *  DA8S(A<7) )  +  DABS(A<8>) 

IF  <D4.EQ.0>00  TO  15 
C. .  .GET  COEFFICIENTS  OF  LINEAR  EQN  BY  ELIMINATING  NONLINEAR  TERM 
B< 1 )  =  A  <  6 ) *A  <  7 )  -  A(5)*A(S) 

B( 2)  =  A(4)*A<7)  -  A(3)*A(S> 

B(3)  =  A  <  2 ) *A  <  7 )  -  A(l> *A ( 3 > 

IF(B(1). EQ. 0. 0  .AND.  B(2) . EQ. 0. 0)G0  TO  20 
M»1 

BIGB  »  DABS ( B  <  1  > ) 

IF<  BIGB. LT. DABS<  B<2) )M  =,  2 
D4  a  B(M) 

DO  1  I  a  1,3 

1  B ( I )  =  B<I)/D4 

L  =  7 

BIGA  a  DABS(A(7) ) 

IF( BIGA. LE. BABS< A< 8 ) )L  *  8 
D4  a  A(L) 

DO  2  I  =  1,3 

2  A(I)  =  A(I)/D4 

IF< L. EQ. 7 )  A<2)  a  A<1) 

IF(L. EQ. 7)  A(4)  a  a<3) 

IF(L.EQ.7)  A(6)  =  A<5) 

IF (M. EQ. 2)  GOTO  5 
IF  ( B  ( 2 ) . EQ . 0 . 0 ) GOTO  10 
BB  =  0(3)  +  B < 2 )  *A  ( 6 )  -  A(4))/B(2) 

CC  =  (B(3)*A(6)  -  A ( 2 >  > /B ( 2 ) 

DISC  =  BB*BB  -  4. *CC 
IF (DISC. LT. 0. 0)  RETURN 
X2  =  ( -BB  +  DSQRT (DISC ) ) /2. 

X4  =  -X2  -  BB 
XI  =  -B ( 2 ) *X2  -  B ( 3 ) 

X3  =  -B ( 2 ) *X4  -  B ( 3 ) 

RETURN 

10  DIV  =  A<4 )  -  B ( 3 ) 

IF ( DI V. EQ. 0)00  TO  20 
C... THERE  IS  ONLY  ONE  ROOT 
XI  a  -  BO) 

X2  =  ( A ( 6  >  *B ( 3 )  -  A ( 2 ) ) / D I V 
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15 

20 


RETURN 

IF<BU>. EG). 0.0)00  TO  ll 

BB  *  <B(3)  +  B< 1 )*A<4 )  -  A(6))/B(l) 

CC  *  <B(3)*A(4)  -  A(2))/B<1) 

DISC  =  BB*BB  -  4.*CC 

IF (DISC.LT. 0.0) RETURN 

Xi  =  <-BB  +  DSQRT<DISO )/2. 

X3  =  -XI  -  BB 

X2  =  — B< 1 )*X1  ~  B<3) 

X4  a  -B  ( 1 ) *X3  -  B<3> 

RETURN 

DIV  *  B<3)  -  A<6> 

IF<DIV.EQ.O)GOTO  20 
. ONLY  ONE  ROOT 
X2  ■  -B  <  3 ) 

XI  *  <A<2)  -  A<4)*8(3> )/DIV 

RETURN 

CONTINUE 

00  TO  111 

CONTINUE 

RETURN 

END 
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PACE  5 


MJ J  .  112. M21 , M22  = 

.  347700D+05— . 7303000+05- . 730300D+05  . 3503000+07 
LAST  6  ELEMENT?  OF  DIAG.  DAMPING  MATRIX  = 

.1414000-01  .1329000-01  . 284600D-01  .3331000-01  .3793000-01  .9613000-01 
LAST  3  ELEMENS  OF  DIAG.  STIFF.  MATRIX  = 

.2190000-01  .3290000-01  .2000000+01  .3345000+01  .3100000+01 
.1463000+02  .7741 000+02  .9241000+02 
PH  II  = 

. i 000000+0 1  . 1 000000+0 1  - .  91 2000D-03- . 1 300000-03- . 1 40000D-03 
. 4070000-03  . 2350000-02  . 7 1 30000-02- . 3000000-07- . 1 720000-04 
PH  1 2  = 

.1000000+01  . 100000D+01  . 344000D— 02  .5430000-04  .3600000-02 
. 3470000-02  . 3950000— 02- . 246000D-02-. 1 15000D-03-. 4130000-06 


CURRENT  RUN  HAS  F  =  0 


INPUT  THE  ROW  AND  COLUMN  OF  ALPHA 

1  ,  1 


INPUT  THE  ROW  AND  COLUMN  OF  BEfA 


1.2 


ENTER  VALUES  OF  SIG  AND  OMEGA 

0.  .  .  1 


ENTER  VALUES  OF  DELTA  SIG  AND  DELTA  OMEGA 

0.  .  .  1 


ENTER  VALUES  OF  FINAL  SIG  AND  FINAL  OMEGA 
1 ...  25 

SIG. W. ALPHA, BETA  =  .0Q00000D+01  . 1000000D+00  . 1534355D-02  -.34237320+03 

X(1,1),Y(1,1),X(2,2).Y<2»2)  ...  * 

. 11D+02-. 150-05  . 69D+01-. 12D-03  . 33D-03  . 1SD-07  . 73D-03  .320-26  . 20D-01  . 14D-04 
.330-01  . 13D-04  . SID-01  . 23D-04  .350-03  . 92D-07  .770+00  . 3S0-04-. 62D-09  . 14D-12 
SIG. W. ALPHA, BETA  =  .00000000+01  .20000000+00  .2109157D-01  -. 12677S4D+04 
X< 1. 1),Y(1, 1),X(2,2).YC2,2)  ...  = 

.  420+02- .  42D-04-. 5SD+02-. 36D-02-. 16D-02  . 100-06  .430-03  .000+01  .200-01  . 230-04 
. 33D-01  . 37D-04  .310-01  .570-04  .350-03  .130-06  . 77D+00  .130-03  . 74D-09-. 15D-12 
ENTER  1  TO  CONTINUE 
STOP 

:ai  darpa  sample  run 


o  o  o 
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DOUBLE  PRECISION  MASSC  4) .  C( 6) . STIFF(S) . PHI 1 ( 10) . CC( 10) 

DOUBLE  PRECISION  SIG. W, DSIG. DW, FSIG.FW, AR, AI , BR, BI , CR. Cl . DR. DI 
DOUBLE  PRECISION  ALPHA, BETA, DELR, DELI . DELT. S2R.S2I . T, SUMR, SUM I 
READ <6, 100 )C 
READ  ( 6 , 10OSTIFF 
READ (6, 100 ) PHI 1 
READ (6, 100JCC 
READ (6, 100) MASS 
100  FORMAT ( 10D1 1.4) 

WRITE<  5. 200 ) MASS 

200  FORMAT < '  Ml  1 , M12, M21 . M22  ■  ',/,4D12.6) 

CC( 1)  *  O.ODO 

CC  <  2 )  =  O.ODO 
DO  16  I  «  3,3 
K  *  1-2 
CC(I)  *  C<K) 

16  CONTINUE 

WRITE<5, 201 ) 

201  FORMAT < '  LAST  S  ELEMENTS  OF  DIAG.  DAMPING  MATRIX  »  ') 

WRITE <5, 202) <CC< I), 1-1.3) 

202  FORMAT <5D 15.  7) 

WRITE <5, 203) 

203  FORMAT ("  LAST  S  ELEMENTS  OF  DIAG.  STIFF.  MATRIX  «  ') 

WRITE (5, 202) STIFF 

WRITE<5»  204) 

204  FORMAT < "  PHI1  -  ') 

WR I TE ( 5 , 202 ) PH I 1 

1  WRITE<5,  300) 

300  FORMAT <//  '  ENTER  VALUES  OF  SIG  AND  OMEGA  '/) 

READ <  4, 301 )  SIG, W 

301  FORMAT  <  2D 16.3) 

WRITE (5, 302) 

READ (4, 301 )DSIG,DW 
WRITE(5, 303) 

302  FORMAT (//  '  ENTER  VALUES  OF  DELTA  SIG  AND  DELTA  OMEGA  '/) 

303  FORMAT <//  '  ENTER  VALUES  OF  FINAL  SIG  AND  FINAL  OMEGA  '/) 

READ <4, 301 )FSIG.FW 

2  AR  =  l.ODO 

AI  =  O.ODO 

DELT  “MASS ( 1 ) *MASS ( 4 >  -  MASS(2)*MASS(3) 

...TRYING  TO  PREVENT  OVERFLOW 

DELR  =  l.ODO 
DELI  =  O.ODO 

CALL  CMULT (SIG, W, SIG, W, S2R , S2 I ) 

CALL  CD I V  <  S2R , S£ 1 , 1 . ODO , 0 . ODO , AR , A I ) 

T  =  ( MASS < 1 >  -  2 . DO*MASS ( 2 )  +  MASS( 4 > ) /DELT 
SUMR  =  T*AR 
SUM I  =  T*AI 
DO  5  I  a  1 ,  S 

CALL  QUAD ( CC ( I ) , ST I FF ( I > , BR , B I , S I G , W ) 

CALL  CD I V  <  BR , BI , 1 . ODO , 0 . ODO , CR ,  C I ) 

SUMR  =  PH 1 1 ( I +2 ) *PH 1 1 ( I +2 ) *CR  +  SUMR 
SUM I  =  PHI1 ( I+2)*PHI1 ( I+2)*CI  +  SUM I 
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3  CONTINUE 

CR  -  -SUMR 
Cl  »  -SUMI 

CALL  CMULT(CR,CI,SIG,W,BR,BI> 

AR  »  OELR 
AI  »  DELI 

CALL  SOLVE 1 < AR, AI , BR* BI , CR, Cl , ALPHA, BETA) 

DELR  *  BR*CI  -  BI*CR 

WRITE (4, 21 1)SIG,W,DELR, ALPHA, BETA 

WRITE (S, 212 >SIG,W. ALPHA, BETA 

211  FORMAT ( *  SIG,W,DET, ALPHA, BETA  -  ' , 2DS. 2, 3D12. 3) 

212  FORMAT (3016. S) 

SIG  =*  SIG  +  DSIG 
W  =  W  +  DW 

IF  (  (SIG.LE.FSIG). AND. (W.LE.FW)  >  GO  TO  2 
WRITE<4,312) 

312  FORMAT ( '  ENTER  1  TO  CONTINUE'/ ) 

READ (5,213)11 

213  FORMAT ( 12) 

IF( II.EQ. 1 >G0  TO  1 
STOP 
END 
C 
C 

SUBROUTINE  CMULT < X , Y, U, V, AR, AI > 

DOUBLE  PRECISION  X , Y, U, V, AR, AI 
AR  =  X*U  -  Y*V 
AI  *  X*V  +  Y*U 
RETURN 
END 
C 
C 

SUBROUTINE  CDIV( X , Y, U, V, AR, AI ) 

DOUBLE  PRECISION  X, Y, U, V, AR, AI , 0 
D  *  X*X  +  Y*Y 
IF  (D.EQ.O.ODO)  RETURN 
AR  =  <X*U  +  Y*V)/D 
AI  =  <X*V  -  Y*U)/D 
RETURN 
END 
C 

SUBROUT I NE  QUAD ( A , B , AR , A I , S I G , W ) 

DOUBLE  PRECISION  A, B, AR, AI , SIG, W 
AR  =  SIG*SIG  -  W*W  +  A*SIG  +  B 
A I  =  2.  ODC)*W*SIG  +  A*W 
RETURN 
END 
C 


SUBROUT I NE  SOLVE 1 ( AR , A I , BR, B I , CR, C I , X 1 , X2 ) 
DOUBLE  PRECISION  AR, A I , BR, BI , CR, Cl , XI , X2, DET 
DET  =  BR*CI  -  BI*CR 
I F ( DABS ( DET ) . LT . 1 . OD-37 ) GO  TO  1 0 
XI  =  -Cl /DET 
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X2  »  BI/DET 
RETURN 

10  WRITE<5,100> 

100  FORMAT ( '  TROUBLE,  LET  ALPHA  -  BETA  ■O') 
XI  »  0. 000 
X2  ■  O.ODO 
RETURN 
END 


LDARPA  3 


9  7 


ldarpa 


M11.M12.M21.M22  » 

. 347700D+05-. 730300D+03-. 730S00D+05  . S30600D+07 
LAST  3  ELEMENTS  OF  DIAG.  DAMPING  MATRIX  ■ 


.OOOOOOOD+Ol 

.3331000D-01 


.OOOOOOOD+Ol 

.37980000-01 


LAST  3  ELEMENTS  OF  DIAG.  STIFF. 


. 2190000B— 01 
. 1468000D+02 
PHI  1  » 

. 1000000D+01 
-. 4070000D-03 


. 3290000D-01 
. 7741000D+02 

. 1000000D+01 
. 23500000-02 


. 1414000D-01 
. 9613000D-01 
MATRIX  - 
.20000000+01 
.92410000+02 

-.91 20000D-03 
.71300000-02 


1329000D-01 

334S000D+01 

1300000D-03 

30000000-07 


ENTER  VALUES  OF  SIG  AND  OMEGA 

0.  .  .  1 


ENTER  VALUES  OF  DELTA  SIG  AND  DELTA  OMEGA 

0.  ,  .  1 


ENTER  VALUES  OF  FINAL  SIG  AND  FINAL  OMEGA 

1. . 1.01 


S I G , W . DET , ALPHA , BETA 
SIG, W, DET, ALPHA, BETA 
S I G ,  W , DET , ALPHA , BETA 
SIG,  14, DET, ALPHA, BETA 
SIG,  14,  DET,  ALPHA,  BETA 
SIG, W, DET, ALPHA, BETA 
S I G ,  W , DET , ALPHA , BETA 
SIG,  Ul,  DET,  ALPHA,  BETA 
SIG, U, DET, ALPHA, BETA 
SIG, W, DET, ALPHA, BETA 
ENTER  1  TO  CONTINUE 


. 00D+01  . 10D+00 
. 0QD+01  . 20D+00 
.000+01  . 30D+00 
.000+01  . 40D+00 
. 00D+01  . 50D+00 
.000+01  . 60D+00 
. 00D+01  . 70D+00 
. 00D+01  . 80D+00 
. OOD+O 1  . 90D+00 
. 00D+01  . 100+01 


- . 8506 1 D-06 
-. 124430-06 
— . 33160D-07 
-. 142220-07 
— . 70478D— 08 
-. 39263D— 08 
-. 23623D-08 
-. 1 49S4D— 08 
-. 98573D-09 
-. 66433D— 09 


. 13344D-02 
.210920-01 
. 11299D+00 
. 37731D+00 
. 96788D+00 
.212920+01 
. 42353D+0 1 
. 78631D+01 
. 1 3954D+02 
. 24027D+02 


STOP 
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23460000-01 


31000000+01 


1400000D-03 

17200000-04 


-.  342S7D+03 
-.  12673D+04 
-. 29210D+04 
-. 53034D+04 
-.  34228D+04 
-. 123620+05 
-.  17214D+05 
-.23107D+05 
*-. 30216D+05 
-.387920+03 
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.000000000+01 
.000000000+01 
.000000000+01 
.000000000+01 
.000000000+01 
. OOOOOOOOD+O 1 
.000000000+01 
. OOOOOOOOD+O 1 
.000000000+01 
. OOOOOOOOD+O 1 
.000000000+01 
. OOOOOOOOD+O 1 
.000000000+01 
. OOOOOOOOD+O 1 
. 000000000+01 
. 000000000+01 
.000000000+01 
. OOOOOOOOD+O 1 
. OOOOOOOOD+O 1 
. OOOOOOOOD+O 1 
. OOOOOOOOD+O 1 
. OOOOOOOOD+O 1 
. OOOOOOOOD+O 1 
.000000000+01 
. 000000000+01 
. OOOOOOOOD+O 1 
. OOOOOOOOD+O 1 
. OOOOOOOOD+O 1 
. 000000000+01 
. OOOOOOOOD+O 1 
. 000000000+01 
. OOOOOOOOD+O 1 
. OOOOOOOOD+O 1 
. OOOOOOOOD+O 1 
. OOOOOOOOD+O 1 
. OOOOOOOOD+O 1 
. OOOOOOOOD+O 1 
.000000000+01 
. OOOOOOOOD+O 1 
.000000000+01 
. OOOOOOOOD+O 1 
. OOOOOOOOD+O 1 
. OOOOOOOOD+O 1 
. OOOOOOOOD+O 1 
.000000000+01 
. OOOOOOOOD+O 1 
. OOOOOOOOD+O 1 
. OOOOOOOOD+O 1 
. OOOOOOOOD+O 1 
. OOOOOOOOD+O 1 
. OOOOOOOOD+O 1 
. OOOOOOOOD+O 1 
. OOOOOOOOD+O 1 
.  OOOOOOOOD+O 1 
. OOOOOOOOD+O 1 


OMEGA 


. 10000000D-01 
. 20000000D-Q1 
. 30000000D-01 
. 40000000D-01 
. 50000000D-0 1 
. 60000000D-01 
. 70000000D-01 
.300000000-01 
. 90000000D-01 
. 10000000D+00 
. 11000000D+00 
. 12000000D+00 
. 13000000D+00 
. 1 4000000D+00 
. 1 4500000D+00 
. 1 45 1 OOOOD+OO 
. 1 4520000D+00 
. 1 4530000D+00 
. 1 4540000D+00 
. 1 4550000D+00 
. 14560000D+00 
. 1 4570000D+00 
. 1 4580000D+00 
. 1 4590000D+00 
. 1 4600000D+00 
. 1461 OOOOD+OO 
. 1 4620000D+00 
. 1 4630000D+00 
. 1 4640000D+00 
. 1 4650000D+00 
. 1 4660000D+00 
. 1 4670000D+00 
. 1 4680000D+00 
. 1 4690000D+00 
. 1 4700000D+00 
. 1471 OOOOD+OO 
. 1 4720000D+00 
. 14730000D+00 
. 1 4740000D+00 
. 1 4750000D+00 
. 14760000D+00 
. 14770000D+00 
. 14730000D+00 
. 14790000D+00 
. 14300000D+00 
. 15000000D+00 
. 1 6000000D+00 
. 170000000+00 
. 1 3000000D+00 
. 19000000D+00 
. 20000000D+00 
. 2 1 0000000+00 
. 2 2 0 0 OOOOD+OO 
. 23000000D+00 
. 24000000D+00 
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. 1 45648660-06 
. 233254640-05 
. 11327427D-04 
.374636190-04 
. 9 1 77333 1 D-04 
.  19112995D-03 
. 356145350-03 
. 6 1 230523D-03 
. 99133676D— 03 
. 1 53435470-02 
. 230034 10D-02 
. 3394S322D-02 
.511 240520—02 
. 9721S399D— 02 
. 5939121 ID-01 
. 693390 19D-01 
. 33464236D-01 
. 1 029SS32D+00 
. 1 326756 1B+00 
. 131 64666D+00 
.27261 384D+00 
. 47693608D+00 
.  U326957D+01 
. 67774S53D+0 1 
. 1 965578 1D+02 
. 1 2399394D+0 1 
. 37003492D+00 
. 1 630897 1D+00 
. 85934 104D-01 
. 50050446D— 0 1 
. 31009480D— 01 
. 1 998945 ID-01 
. 13213922D— 01 
. 38632056D-02 
. 59S04427D-02 
. 40275020D-02 
. 26351706D-02 
. 1 75572 17D-02 
. 1 1 122133D— 02 
. 67068735D— 03 
. 374 1 4568D-03 
. 133073130-03 
. 6954 1 03SD-04 
. 134SS117D-04 
. 29S5S543D-06 
. 1 7S00S4SD-02 
. 67900350D-02 
. 99206569D— 02 
. l 3 1 69S96D— 0 1 • 
. 168482960-01 
.210915740-01 
. 259980 17D-01 
.31  o5c-863D— 0 1 
. 381573470-01 
. 455933500-01 


F( 1*2) 


-.334353190+01 
-. 133801350+02 
-.30123155D+02 
-.536203650+02 
- . 339 1 1 7S6D+02 
-. 121 08379D+03 
-. 1 6526630D+03 
-.216668540+03 
-.275648610+03 
-.342873220+03 
-.41 974390D+03 
-. 50980860D+03 
- . 62547862D+03 
-.36231 907D+03 
-.214002730+04 
-.231 0924 1D+04 
-.252630440+04 
-.230626260+04 
-. 31S51415D+04 
-.372633300+04 
— . 45656735D+04 
-. 60339280D+04 
-. 93064754D+04 
-. 2276469 1D+05 
. 33767323D+05 
. 97370053D+04 
. 53 191 835D+04 
.353131250+04 
. 25633299D+04 
. 1 9562509D+04 
. 1 5398083D+04 
. 1 2362853 D+04 
. 1 005 1 549D+04 
.323211330+03 
. 6762098 1D+03 
. 55492 147D+03 
.453104560+03 
. 366386250+03 
. 29161 187D+03 
. 22644330D+03 
. 1 69 1 3295D+03 
. 1 13309350+03 
. 729 1 6577D+02 
.321 12937D+02 
-. 477790500+01 
-.368892750+03 
-. 72027024D+03 
-.370358620+03 
100249310+04 
-. 1 1335024D+04 
-. 1 2677S3SD+04 
-. 140701610+04 
1 5520031 D+04 
-. 17C32190D+04 
1 3C' 08 998 D+04 


.*» 


PRODUCTION  RUN  FOR  SIGMA 


0 


PAGE  1 


SIGMA 


. OOOOOOOOD+Ol 
. OOOOOOOOD+Ol 
. OOOOOOOOD+O 1 
. OOOOOOOOD+O 1 
. OOOOOOOOD+O 1 
. OOOOOOOOD+O 1 
. OOOOOOOOD+O 1 
. OOOOOOOOD+O 1 
. OOOOOOOOD+O 1 
. OOOOOOOOD+O 1 
. OOOOOOOOD+O 1 
. OOOOOOOOD+O 1 
. OOOOOOOOD+O 1 
. OOOOOOOOD+O 1 
. OOOOOOOOD+O 1 
. OOOOOOOOD+O 1 
. OOOOOOOOD+O 1 
. OOOOOOOOD+O l 
. OOOOOOOOD+O 1 
. OOOOOOOOD+O 1 
. OOOOOOOOD+O 1 
. OOOOOOOOD+O 1 
. OOOOOOOOD+O 1 
. OOOOOOOOD+O 1 
. OOOOOOOOD+O 1 
. OOOOOOOOD+O 1 
. OOOOOOOOD+O 1 
. OOOOOOOOD+O 1 
. OOOOOOOOD+O 1 
. OOOOOOOOD+O 1 
. OOOOOOOOD+O 1 
. OOOOOOOOD+O 1 
. OOOOOOOOD+O 1 
. OOOOOOOOD+O 1 
. OOOOOOOOD+O 1 
. OOOOOOOOD+O 1 
.OOOOOOOOD+Ol 
. OOOOOOOOD+O 1 
. OOOOOOOOD+Ol 
.OOOOOOOOD+Ol 
.OOOOOOOOD+Ol 
. OOOOOOOOD+Ol 
•OOOOOOOOD+Ol 
•OOOOOOOOD+Ol 
.OOOOOOOOD+Ol 
.OOOOOOOOD+Ol 
•OOOOOOOOD+Ol 
. OOOOOOOOD+Ol 
. OOOOOOOOD+Ol 
. OOOOOOOOD+Ol 
. OOOOOOOOD+Ol 
.  OOOOOOOOD+O'l 
. OOOOOOOOD+O 1 
. OOOOOOOOD+Ol 
. OOOOOOOOD+O 1 


OMEGA 


. 25000000D+00 
. 26000000D+00 
. 27000000D+00 
. 28000000D+00 
. 287 OOOOOD+OO 
. 287 1 OOOOD+OO 
. 23720000D+00 
. 23730000D+00 
.287400000+00 
. 28750000D+00 
. 28760000D+00 
. 28770000D+00 
. 28780000D+00 
. 28790000D+00 
. 28800000D+00 
.2881 OOOOD+OO 
. 28820000D+00 
. 28830000D+00 
. 28840000D+00 
. 30000000D+00 
.SOOOOOOOD+OO 
. 70000000D+00 
. 90000000D+00 
. 11000000D+01 
. 13000000D+01 
. 1 3000000D+0 1 
. 1 7000000D+0 1 
. 1 9000000D+0 1 
.209000000+01 
.209100000+01 
.209200000+01 
. 20930000D+0 1 
. 20940000D+0 1 
. 20950000D+0 1 
. 20960000D+0 1 
. 20970000D+01 
. 20980000D+0 1 
. 20990000D+0 1 
.220000000+01 
. 22500000D+0 1 
. 23000000D+0 1 
. 23500000D+01 
. 24000000D+01 
. 24500000D+01 
.250000000+01 
. 25500000D+0 1 
. 26000000D+01 
. 26500000D+01 
. 27000000D+01 
. 27500000D+0 1 
.280000000+01 
. 28500000D+01 
. 29000000D+0 1 
. 29500000D+0 1 
. 300000000+01 
95 


F(l.l) 


.540686400-01 
. 637096 16D-01 
. 74718126D-01 
. 87826588D-01 
.114079500+00 
. 1 1689545D+00 
. 12058962D+00 
. 12567 636D+00 
. 133166580+00 
. 145351720+00 
. 168730420+00 
.231258580+00 
.751199120+00 
. 17327799D-01 
. 23279390D-01 
. 46085 1 80D-0 1 
.585562850-01 
. 662S7989D-01 
.714805910-01 
.112993330+00 
. 96788292D+00 
. 42353357D+0 1 
. 1395381 1D+02 
. 4080905 1D+02 
. 1 2986081 D+03 
. 4067 1 942D+03 
. 1 9663688D+04 
. 97089372D+04 
. 488340 18D+07 
. 5407 4079D+07 
. 57967006D+07 
. 59705091D+07 
. 58879446D+07 
. 55700729D+07 
. 508803 17D+07  . 
. 45282139D+07 
. 39620073D+07 
. 34341 702D+07 
. 27463862D+05 
. 1 377 4664D+05 
-.  35686 149D+04 
. 60248 148D+04 
. 45903 148D+04 
. 37084945D+04 
. 31377357D+04 
. 27632936D+04 
. 252764 14D+04 
. 24052460D+04 
. 23941849D+04 
. 25192944D+04 
. 285 1 9092D+04 
. 3573031 5D+04 
.51 9 5 o 3 3 o D + 0 4 
.  9604267 4D-r 04 
. 23432593D+05 


F(l,2) 


-. 20256965D+04 
-.  21978706D+04 
-.237904400+04 
-.257800780+04 
-. 29370909D+04 
-.297310400+04 
-.301970120+04 
-.308271610+04 
-. 3 1 7323380+04 
-.331521920+04 
-.357137760+04 
-.418163880+04 
-. 75365544D+04 
.114462890+04 
-.132671160+04 
-. 1 8666766 D+04 
—.210413160+04 
-.22382211 D+04 
-.232474660+04 
-.292103370+04 
— . 84228344D+04 
-.  172138910+05 
-.302163750+05 
-.491941790+05 
— . 78063751D+05 
-. 1 2396795D+06 
-.228541650+06 
—.371697290+06 
-.489741320+07 
-.373527320+07 
— .2198S565D+07 
— .. 40837 4590+06 
. 1 4240560D+07 
.307184580+07 
. 43781 987D+07 
• 52927990D+07 
. 58501 557D+07 
.612527950+07 
.  7920 1 906D+06 
.537901550+06 
. 40357793D+06 
. 3 1 928990D+06 
.260481930+06 
.21 625722D+06 
. 1 S099226D+06 
. 15143085D+06 
. 12547728D+06 
. 10161 035D+06 
. 7S5363S2D+05 
. 54902433D+05 
•28919939D+05 
-.23796301 D+04 
— . 45 i 14363D+05 
-. 11541321D+06 
-.  2744S942D+06 


PRODUCTION  RUN  FOR  SIGMA  *  0 


PAGE 


SIGMA 


OMEGA 


F(l.l> 


F( 1.2) 


. OOOOOOOOD+O 1  . 305000000+0 1 
.000000000+01  .310000000+01 
. 000000000+0 1  .31 3000000+0 1 
. 000000000+0 1  . 320000000+0 1 
. 000000000+0 l  . 325000000+0 1 
. 000000000+0 1  . 330000000+0 1 
. OOOOOOOOD+O 1  . 333000000+0 1 
. 000000000+0 1  . 340000000+0 1 
. 000000000+0 1  . 345000000+0 1 
. 000000000+0 1  . 350000000+0 1 
. 000000000+0 1  . 35500000D+0 1 
. 000000000+0 1  . 34000000D+0 1 
. 000000000+0 1  . 345000000+0 1 
. 000000000+0 1  . 37000000D+0 1 
.000000000+01  .375000000+01 
. 000000000+0 1  . 380000000+0 1 
. OOOOOOOOD+O 1  . 385000000+0 1 
. OOOOOOOOD+O 1  . 39000000D+0 1 
. OOOOOOOOD+O 1  . 39500000D+0 1 
. OOOOOOOOD+O 1  . 40000000D+0 1 
. OOOOOOOOD+O 1  . 40500000D+0 1 
. OOOOOOOOD+O 1  .41 OOOOOOD+O 1 
. OOOOOOOOD+O 1  .41 5000000+0 1 
. OOOOOOOOD+O 1  . 42000000D+0 1 
. OOOOOOOOD+O 1  . 42500000D+0 1 
. OOOOOOOOD+O 1  . 43000000D+0 1 
. OOOOOOOOD+O 1  . 43500000D+0 1 
. OOOOOOOOD+O 1  • 440000000+0 1 
.  OOOOOOOOD+O  1  44500000D+0 1 
. OOOOOOOOD+O 1  . 430000000+0 1 
. OOOOOOOOD+O 1  . 45500000D+0 1 
. OOOOOOOOD+O 1  . 440000000+0 l 
. 000000000+0 1  . 44500000D+0 1 
. OOOOOOOOD+O 1  . 470000000+0 1 
. OOOOOOOOD+O 1  . 47500000D+0 1 
. OOOOOOOOD+O 1  . 480000000+0 1 
. OOOOOOOOD+O 1  . 48500000D+0 1 
.-OOOOOOOOD+O  1  . 49000000D+0 1 
. OOOOOOOOD+O 1  . 49500000D+0 1 


.321814490+04  -. 47218240D+04 
. 44 1 342930+03  . 344337390+04 
. i 44574740+05  . 30030844D+04 
.451581 040+04  . 204983490+04 
. 389424790+04  . 1 58289780+04 
. 271141030+04  . 127082580+04 
. 20703 151 D+04  . 1 0430872D+04 
. 14805598D+04  .848429340+03 
. 1 4233 1 800+04  . 72244497D+03 
.124284250+04  .594971740+05 
.111011 35D+04  . 48592249D+03 
. 1 00878990+04  . 3854 1 4230+03 
.929044180+03  . 29283 134D+05 
. 344754450+03  . 2044 1 306D+05 
.811801 440+03  . 1 24792070+05 
. 74742820D+03  . 47014446D+04 
.729491140+03  -.27451430D+04 
. 497 1 875 1 D+03  - . 9978405 1 D+04 
.643882420+03  -. 149839 19D+05 
.  443994900+03  - .  238 1 7444D+05 
.421 93398D+03  - . 303037 45D+03 
. 40222374D+03  - . 37082 1 30D+05 
. 5845 1 395D+03  - . 4355803 1 D+05 
. 54S49553D+03  - . 4995334 1 D+05 
. 55393554D+03  - . 542824440+05 
. 3404390 1 D+03  - . 42538247D+0S 
. 52344484D+03  - . 48790392D+03 
. 51721870D+03  -.749897700+05 
. 50484759D+03  -.811 42973D+03 
. 497235790+03  -. 87317581D+05 
. 48830 1 43D+03  - . 93459833D+05 
. 47997494D+03  - . 99395244D+05 
.4721 95 1 9D+03  -.105728790+04 
. 44490949D+03  -.11 1S4481D+04 
. 45807230D+03  -.11 300724D+04 
.451 44330D+03  -.1241 5949D+04 
. 44553455D+03  - . 1 3032529D+04 
.439870800+03  -. 13450497D+04 
. 434448 1 1 D+03  - . 1 427 0735D+04 
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PRODUCTION  RUN  FOR  3IGMA  *  0 


PAGE 


SIGMA 


. OOOOOOOOD+Ol 
. OOOOOOOGD+O 1 
.000000000*01 
.  000000000*01 
. 000000000*01 
.  000000000*01 
.  000000000*01 
. 000000000*01 
. 000000000*0 1 
. 000000000*01 
. 000000000*01 
.000000000+01 
. 0Q000000D+0 1 
.000000000*01 
. 000000000*01 
. 000000000+0 1 
. 000000000+01 
.000000000+01 
. 000000000+01 
. 000000000+01 
. OOOCOOOOD+O 1 
. 000000000+01 
.000000000+01 
. 000000000+01 
.000000000+01 
. 000000000+01 
.000000000+01 
.000000000*01 
. 000000000+01 
. 000000000+0 1 
.000000000+01 
. OOOOOOOOD+Ol 
.000000000+01 
.000000000+01 
. 000000000+01 
. 000000000+01 
.000000000+01 
. 000000000+01 
. 000000000+01 
. 00000000D+01 
. 000000000+0 1 
. OQOCOOOOD+O 1 
. 000000000+0 1 
. OOOOOOOOD+Ol 
. OOOOOOOOD+Ol 
.  OOOOOOOOD+Ol 
.  OOOOOOOOD+0 1 
. OOOOOOOOD+Ol 
.  OOOOOOOOD+Ol 
.  OOOOOOOOD+Ol 
.  OOOOOOOOD+0 1 
.  OGOOOOCGD+G 1 
.  OOOOOOOOD+Ol 

.  OOOOOOOOD+Ol 


ENLARGEMENT 


OMEGA 


. 1 OOOOOOOD-O 1 
.  1 1000000D-01 
.  12001 ’0000-0 1 
.  1 3000000 D—0 1 
. 14000000D-01 
. ISCOOOOOD-Ol 
. 1 6000000D-0 1 
. 17000000D-01 
. ISOOOOOOD-Ol 
. 190000000-01 
.200000000-01 
.210000000-01 
. 22000000D— 01 
.230000000-01 
. 24000000D— 0 1 
. 2S000000D— 0 1 
.260000000-01 
■  .270000000-01 
. 2S000000D-0 1 
.290000000-01 
. 30000000D-0 1 
. 31000000D-01 
. 32000 OOOD—O 1 
. 330000000-01 
. 34000000D-0 1 
.330000000-01 
. 36000000D— 0 1 
.  370000000-01' 

. 380000000-01 
.390000000-01 
. 40000000D-0 1 
.41 OOOOOOD-O 1 
. 42000000D-01 
. 43000000D-0 1 
. 440C OOOOD-Ol 
.450000000-01 
.  14802000D+00 
.  1 4804000D+00 
. 1 4806000D+00 
.  1480S000D+00 
.  1481 00000+00 
. 1 43 1 20000+00 
. 148140000+00 
. 148160000+00 
. 148180000+00 
. 148200000+00 
. 287924000+00 
. 287924200-00 
. 287924400+00 
. 287924600+00 
. 287924800— 00 
. 287925 0 0 0-00 
. 237925200+00 
.287923400+00 
. 287925600-00 
98 

JF  STAc-Lc.  REG  I  Of. 


F  ( 1 .  1  ) 


. 1 4564866D-06 
.213237920-06 
. 302037830-06 
.41 6073290—06 
. 339689040-06 
. 737630070-06 
. 934982040— 06 
. 121718180-03 
.  133002000-05 
.  189964160-03 
. 233234640-03 
.  .*86556990—05 
. 341398400-05 
.408129760-03 
. 483945770-03 
. 56987 49SD— 05 
. 666781 94D— 03 
. 773567260-05 
.8971 677 ID-05 
. 1 0323564D-04 
.118274270-04 
. 1 3487727D— 04 
. 153172910-04 
. 1 73273 12D-04 
. 193293560-04 
.219333560-04 
. 245576200-04 
.274088310-04 
.  3050 204 7D-04 
. 33S50704D-04 
. 374686 19D-04 
.41 36999 ID-04 
. 435694040-04 
. 30081 S29D— 04 
. 54922623D-04 
.601075550-04 
. 1S017429D-05 
.450817250-03 
. 835050050-05 
. 132651350-04 
. 191920320-04 
. 2c0 7 4*  780—04 
.  33:3588840—04 
. 424946000-04 
.519337360-04 
.621309970-04 
. 242043910-05 
. 542763460-03 
. 9603377 i 0-05 . 
. 1 49465450— 04 
. 214232810-04 
•  2:‘06 1 9640— '.**■ 


. 583826860—04 


(SIGMA  =  0) 


F< 1,2) 


-. 334353190+01 
-.404585670+01 
-. 48 1506390+01 
-.565121700+01 
-. 655433100+01 
-. 752442190+01 
-.336150710+01 
-. 966560520+01 
-. 108367360+02 
— . 12G74922D+02 
- . 1 3380 1 830+02 
-. 147523500+02 
-. 161920410+02 
-.  176986370+02 
-.  192723170+02 
-.209135600+02 
—  •  8*62 1 8480+0*. 
-.243974160+02 
-. 262402980+02 
-.  281505310+02 
-.301281530+02 
-.  321732100+02 
-.342857390+02 
-. 36465737D+02 
-.387134010+02 
-.  41023629D+02 
— .  4341 1523D+02 
— . 45862 1 860+02 
-.4S3S0525D+02 
- . 50966743D+02 
-.53620S65D+02 
-. 5634294 1D+02 
-.591 33042D+02 
-.61991.2390+02 
- . 649 1 7602D+02 
-.679122090+02 
-.  11 73679 1D+02 
-. 1 3565344D+02 
-. 25267295D+02 
-. 31346236D+02 
-. 3S305622D+02 
-.  446*87320+02 
— . 5037S? 2 1 0+02 
— . 569991230+02 
— . 630 1233 10+02 
— . 6892 1 5500+02 
— .  1 656-Z' .  coO+O* 
-.202510560+02 
—  •  8Z‘954z-4*D  +  V».- 
— . 336 173030*02 
— . 40247 1 7 30+02 
- .  -i  68  4  *r  2 1 30+02 


— . 664407 4OD+02 


RAGE  1 


•5  T 


. 000000000+01 
.00000000D+01 


ENLARGEMENT  OF 


OMEGA  F( 1 . 1 > 

. 237923800+00  . 70303230D-04 

. 237926000+00  . 832640730-04 


99 


-  STABLE  ."'Ey I UN  (SlUMA  — 


F(l,2> 

-.729083430+02 

793433070+02 


0) 


PAGE 


SIGMA 


OMEGA 


FU.l) 


F<1,2) 


-.141 400000-0 1  . 1 00000000-02 

-.141 400000-0 1  . 200000000-02 

-.141 400000-0 1  . 300000000-02 

-.141 400000-0 1  . 400000000-02 

-.141 400000-0 i  . 300000000-02 

-.141 400000-0 1  . 600000000-02 

-.141 400000-0 1  . 700000000-02 

-.141 400000-0 1  . 300000000-02 

-.141 400000-0 1  . 900000000-02 

-.141 400000-0 1  . 1 00000000-0 1 

-.141 400000-0 1  . 1 00000000-0 1 

-. 141400000-01  .200000000-01 

-.141 400000-0 1  . 300000000-0 1 

-.141 40C00D-0 1  . 400000000-0 1 

-.141 400000-0 1  . 500000000-0 1 

-.141 400000-0 1  . 600000000-0 1 

-.141 400000-0 1  . 700000000-0 1 

-.141 400000-0 1  . 300000000-0 1 

-.141 40000D-0 1  . 900000000-0 1 

-.141 400000-0 1  . 1 00000000+00 

-. 14140000D— 01  . 110000000+00 

-.141400000-01  .120000000+00 

-.141400000-01  .130000000+00 

-.141400000-01  . 14000000D+00 

-.141 400000-0 1  . 1 30000000+00 

-.141400000-01  .160000000+00 

-.141 400000-0 1  . 1 70000000+00 

-.141 400000-0 1  . 1 30000000+00 

-. 141400000-01  . 190000000+00 

-.141400000-01  . 20000000D+00 

-.141 40000D-0 1  .21 0000000+00 

-.141 400000-0 1  . 220000000+00 

-.141 400000-0 1  . 230000000+00 

-.141 40000D— 0 1  . 24000000D+00 

-.141 400000-0 1  . 250000000+00 

-.141 400000-0 1  . 260000000+00 

-.141400000-01  .270000000+00 

-. 141400000-01  .230000000+00 

-.141400000-01  .290000000+00 

-.141 400000-0 1  . 300000000+00 

-.141400000-01  .310000000+00 

-.141 400000-0 1  . 320000000+00 

-.141 40000D-0 1  . 330000000+00 

-.141 400000-0 1  . 340000000+00 

-.141 400000-0 1  . 330000000+00 

-.141 400000-0 1  . 360000000+00 

-.141 400000-0 1  . 370000000+00 

-.141 400000-0 1  . 330000000+00 

-.141 40000D-0 1  . 390000000+00 

-.141 400000-0 1  . 400000000+00 

-.141 400000-0 1  .41 0000000+00 

-.141 400000-0 1  . 420000000+00 

-.141 400000-0 1  . 430000000+00 

-.141 400000-0 ..  .  440000000+00 

-.141 400000-0 l  . 450000000+00 

100 


-.944393130+03  -. 671194470+01 

- . 944903960+03  -.63121 3000+0 1 

- . 9449 1 3930+03  - . 697924030+0 1 

-.944937100+03  -. 721312900+01 

-.944960490+03  -. 731334910+01 

944939 1 1 0+03  - . 733 140360+0 1 

- . 943022930+03  - . S3 1 530430+0 1 

- . 943062 1 20+03  - . 33 1 703270+0 1 

- .  945 1 06530+03  - . 9333 1 5320+0 1 

-.9431 36370+03  - . 1 0020 1 300+02 

- . 943 1 56370+03  - . 1 0020 1 300+02 

- . 945953030+03  - . 200503290+02 

- . 947362 1 00+03  - . 367359790+02 

-. 949434300+03  -. 602574070+02 

- . 9523 1921 0+03  - . 905 1 234 1 0+02 

- . 936737030+03  - . 1 27622770+03 

- . 962323930+03  -.171 696630+03 

-.971 564600+03  - . 222397950+03 

- . 934735450+03  -.231 43 1 440+03 

-. 100377330+04  -.347349190+03 

-. 104222900+04  -.422616310+03 

■•111 232 1 90+04  - . 506436390+03 

■.12611 4970+04  - . 397699490+03 

•. 133031220+04  -.673733490+03 

•.1541 93 1 00+04  - .  7034239 10+03 

•.124110140+04  -.731949620+03 

.107331 950+04  - . 394050320+03 

. 101023640+04  -. 10161699D+04 

. 973977000+03  -.11 4353700+04 

. 963279270+03  - . 1 27623090+04 

. 93499 1 200+03  -.141 469440+04 

- . 930639 1 60+03  - . 1 339 1 9370+04 

— . 94S340610+03  -.171007320+04 

-. 943323110+03  -. 136734330+04 

- . 950656690+03  -.20313171 0+04 

- . 953242240+03  - . 220303230+04 

-.965413490+03  -. 233127250+04 

-. 933939710+03  -.256420920+04 

- . 996650660+03  - . 2746 1 3390+04 

- . 975355370+03  - . 293300060+04 

-.  962342330+03  -.314279540+04 

- . 957390340+03  - . 33554 1 1 1 0+04 

- . 956006040+03  - . 3574335 1 0+04 

- .  956 1 63500+03  - . 330 1 23920+04 

- . 957 1 03950+03  - . 403459350+04 

-. 953432450+03  -. 427505560+04 

■ . 960 1 47350+03  - . 45227 1 520+04 

. 9620 1 7300+03  - . 477764990+04 

. 964046360+03  - . 303993030+04 

. 966209 1 60+03  - . 530962430+04 

. 963437650+03  - . 553679690+04 

. 970371730+03  -. 537131170+04 

. 973354600+03  -.61 6333400+04 

. 97593 1 720+03  - . 646332930+04 

. 973600 1 90+03  - . 677 1 564 1 0+04 


PARTIAL  PRODUCTION  RUN  FOR  SIGMA  =  -.01414  PAGE  1 


SIGMA 


OMEGA 


F<  1 . 1 ) 


F<1,2> 


-.141 40000D-0 1  . 460000000+00 

-.141 40000D-0 1  . 470000000+00 

-.141 40000D-0 1  . 430000000+00 
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-.141400000-01  .870000000+00  - 
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-.141400000-01  .390000000+00  - 

-.141 400000-0 1  . 900000000+00  - 

-.141 400000-0 1  .91 0000000+00  - 

-.141400000-01  .920000000+00  - 

-.141400000-01  .930000000+00  - 

-.141400000-01  .940000000+00 

-.141400000-01  .950000000+00  - 

-.141 400000-0 1  . 960000000+00  - 

-.141400000-01  .970000000+00  - 

-.141400000-01  .930000000+00  - 

-.141400000-01  .990000000+00  - 

-.141400000-01  .100000000+01 
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141400000-01 
-. 141400000-01 
14140000D-01 
141400000-01 
-. 141 40000D— 0 1 
141400000—01 
14140000D— 01 
14140000D-01 
14140000D-01 
141400000—01 
-. 14140000D-01 
141400000-01 
141400000-01 
141400000-01 
-. 141400000-01 
-. 141400000-01 
-. 141 40000D— 0 1 
-.  141400000-01 
-. 141400000-01 
-.111400000—01 
-. 141400000-01 
-.  141 40000D-0 1 
-. 141 40000D-0 1 
-. 1 41 40000D— 0 1 
-. 141400000-01 
-.  141400000-01 
-. 141 40000D— 0 1 
14140000D— 01 
-. 141400000-01 
141400000-01 
-. 141400000-01 
-. 141400000-01 
-. 141 40000D— 0 1 
-.  141400000-01 
-. 141 40000D— 0 1 
-. 141400000-01 
-. 141400000-01 
-. 141400000-01 
-..141400000-01 
-. 141400000-01 
-. 14140000D— 01 
-. 141400000-01 
-. 141400000-01 
-. 141400000-01 
-. 141400000-01 
-. 141400000-01 
-. 141400000-01 
-. 141400000-01 
-. 141400000-01 
-. 141400000-01 
-. 141400000-01 
-. 141400000-01 
-. 141400000-01 
-. 141400000-01 
-. 141400000-01 


. 101000000+01 
. 10200000D+01 
. 103000000+01 
. 104000000+01 
. 105000000+01 
. 1 0600000D+0 1 
. 1 0700000 D+01 
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. 109000000+01 
. 110000000+01 
. 111000000+01 
. 112000000+01 
.113000000+01 
. 114000000+01 
. 115000000+01 
. 116000000+01 
. 1 1700000D+01 
. 113000000+01 
. 119000000+01 
. 120000000+01 
. 121000000+01 
. 122000000+01 
. 123000000+01 
. 124000000+01 
. 125000000+01 
. 126000000+01 
. 127000000+01 
. 123000000+01 ' 
. 129000000+01 
. 130000000+01 
. 131000000+01 
. 132000000+01 
. 133000000+01 
. 134000000+01 
. 135000000+01 
. 136000000+01 
. 137000000+01 
. 1 38000000+0 1 
. 139000000+01 
. 140000000+01 
. 141000000+01 
. 142000000+01 
. 143000000+01 
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.  146000000+01 
. 147000000+01 
. 1 43000000+0 1 
. 149000000+01 
. 150000000+01 
. 151000000+01 
. 152000000+01 
. 153000000+01 
. 154000000+01 
. 155000000+01 
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-. 131853410+04 
-. 133000820+04 
-. 134179990+04 
-. 135392110+04 
136633460+04 
-. 137920360+04 
-. 139239240+04 
-. 140596600+04 
-. 141994060+04 
-. 143433340+04 
-. 144916270+04 
-. 146444S3D+04 
-. 148021150+04 
-. 149647540+04 
-. 151326490+04 
153060730+04 
-. 154353270+04 
-. 156707420+04 
-. 158626360+04 
-. 160615740+04 
-. 162673770+04 
-. 164321330+04 
-.167049720+04 
-. 169371360+04 
-. 171795130+04 
-. 174332160+04 
-. 176996090+04 
-. 179304770+04 
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-. 139335530+04 
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225224330+04 
-.241323540+04 
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— .  3 2 8 7 6 7 98 D  i-  0  4 
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- . 3030 1 8630+04 
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-. 407143900+05 
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-. 447360930+05 
-.458552920+05 
-. 469460960+05 
-.  480591000+05 
-.  491949230+05 
-.503542130+05 
-.515376670+05 
-.527459390+05 
-.539799400+05 
-.  552403 17D+05 
— . 56527962D+05 
-.578437640+05 
-.591386680+05 
- . 605636770+05 
-.619693620+05 
- . 634033660+05 
-.643304220+05 
-.663373600+05 
-.679306260+05 
-.695113060+05 
-.711326530+05 
-.727951300+05 
— . 7450146SD+05 
-.762542490+05 
-.730565260+05 
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-. 1 17561250+06 
— . 120713310+06 
-. 123940670+06 
-. 1272591 10+06 
-. 130631630+06 
-. 134219920+06 
-. 137335270+06 
-.141 639 1 00+06 
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5.3.  EXTENSION  TO  DIGITAL  SYSTEMS 

The  synthesis  and  analysis  of  the  digital  parameter  plane 
is  very  similar  in  form  to  that  of  the  continuous  technique.1 
In  order  to  understand  It  we  will  consider  the  third  order 
system  that  Is  depicted  in  Fig.  5.6.  The  system  can  be  nor¬ 
malized  as  indicated  so  that  the  results  obtained  will  have 
a  more  general  nature. 

We  are  Interested  in  the  location  of  the  roots  relative 
to  the  unit  circle  (see  FI?.  5.7.).  Roots  lying  Inside  of  the 
unit  circle  are  stable,  roots  lying  outside  of  the  unit  circle 
are  unstable,  and  the  unit  circle  then  represents  one  of  the 
stability  boundaries.  The  other  two  boundaries,  while  actu¬ 
ally  lying  on  the  unit  circle  have  unique  characteristics, 
are  define  at  z  a  ±1. 

We  can  now  make  the  necessary  substitution  for  z  as  we  did 
for  s  and  reduce  the  characteristic  equation  to  two  equations 
in  two  unknowns.  The  form  then  of  the  solution  Is  basically 
the  same  whether  one  is  tal’klng  about  continuous  domain  or 
sampled  domain  systems  (see  Fig.  5.8.).  Also  In  the  digital 
systems  one  can  make  use  of  matrix  formulations  (Sec.  5.2.)  and 
be  able  to  accommodate  high  order  systems  with  relative  ease. 

An  example  of  the  digital  domain  parameter  olane  technique  is 
contained  in  the  following  charts. 
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Fig.  5.8.  COMPLEX  CONJUGATE  BOUNDARY  AT  z  *  e16 

zN  .  eiN0 
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5.4  EXAMPLE 

In  order  to  portray  the  power  of  that  parameter  plane 
technique  we  have  chosen  the  model  as  depicted  on  Fig.  5.9. 

Here  we  have  a  general  digital  control  system  (PD)  as  applied 
to  the  end  of  a  free-free  structure.  The  controller  utilizes 
a  rate  signal  in  the  feedback  loop  that  Is  digitally  Integrated 
to  yield  an  estimate  of  the  vehicle  attitude.  The  parameter 
plane  for  this  configuration  is  given  on  the  following  pages  for 
one  and  then  three  bending  modes  In  the  model. 

The  parameter  plane  for  one  bending  mode  Is  presented  In 
Fig.  5.10.  The  stability  boundaries  and  the  constant  damping 
contours  are  plotted  as  function  of  the  two  control  gains, 
position  (aQ)  and  rate  (aj).  The  stable  region  Is  Indicated 
by  the  cross-hatched  region.  Also  presented  are  contours  of 
constant  closed  loop  damping.  The  designer  can  then  select 
the  values  of  the  two  control  gains  to  achieve  the  stability  and 
response  that  he  desires.  Also,  one  would  be  well  advised  to 
select  a  set  of  gains  that  are  comfortably  within  the  stability 
region  (robustness).  The  plot  on  Fig.  5.11  is  a  "blow-up"  of 
the  origin  of  the  stability  region  to  illustrate  the  movement  of 
the  constant  damping  contours  i-n  the  vicinity  of  the  origin.  It 
is  important  that  the  control  system  analyst  carefully  investigate 
all  regions  of  the  parameter  Diane  to  insure  that  no  "funnies" 
go  unnoticed.  This  fact  will  become  more  apparent  as  we  invest¬ 
igate  the  stability  boundaries  and  contours  for  the  same  control 


configuration  with  the  addition  of  the  next  two  bending  modes. 

This  information  is  contained  on  Fig.  5.12  -  5.13. 

The  parameter  plane  Including  three  bending  modes  Is 

presented  on  Fig.  5.12.  The  scale  utilized  Is  the  same  scale 

as  for  the  one  bending  mode  case  to  emphasize  the  reduction 
In  the  stability  region.  This  graph  forecefully  depicts  the 

requirement  for  adequate  bending  suporesslon  In  high  order  systems. 

The  plot  on  Fig.  5.13  Is  a  more  reasonable  scaling  of  the  stable 

region  for  the  three  bending  mode  case.  The  curves  In  this 

figure  can  be  used  as  guides  in  selection  of  the  appropriate 

gains  to  use  in  not  only  Insuring  stability  but  also  in  providing 

the  requisite  response  characterl s tl cs . 

Again  a  "blow-up"  of  the  region  In  the  vicinity  of  the 

origin  Is  shown  In  Fig.  5.14.  The  motion  of  the  constant  damping 

contours  becomes  more  Interesting  as  the  order  of  the  system 

Is  increased.  Also  we  must  be  very  careful  of  the  increment 

that  Is  used  to  vary  w  t  along  the  constant  damping  contour. 
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ROTE  GAIN 


RATE  GRIM 


wb  1 .00,2.76, 5»40  rad/sec 
T  O.I1626  seconds 
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5.5  SUMMARY 

In  summary,  we  have  presented  a  technique,  valid  for  both 
continuous  and  digital  control  system  design,  that  Is  both 
powerful  In  concept  and  application  but  very  simple  to  under¬ 
stand  and  to  use.  We  have  applied  this  theory  to  a  number  of 
vehicles  and  configurations  over  the  years  and  have  found  It 
to  be  extremely  effective.  The  recent  advancements  under  this 
contract  to  extend  this  technique  to  include  high  order  systems 
hasopenedthe  door  to  exciting  and  Interesting  future  applications 
in  the  analysis  and  design  of  future  large  space  structures. 
However,  further  effort  Is  Indicated  to  handle  cases  whose  char¬ 
acteristic  equations  have  closely-spaced  roots. 
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6. 'I  CONCLUSIONS 

This  report  has  presented  an  overview  of  the  activities  of 
the  Control  Dynamics  Company  In  performing  the  requirements 
of  the  subject  contract.  We  have  presented  a  study  of  active 
structural  control  technology  as  applied  to  Large  Space  Struc¬ 
tures  (LSS).  In  particular,  we  have  presented  the  development 
of  analytical  tools  for  analyzing  the  behavior  of  LSS's  and 
presented  the  development  of  means  to  ensure  LSS  stability  and 
performance.  In  addition  we  have  supported  the  overall  ACOSS 
activities  to  develop  and  understand  a  generic,  unified,  struc¬ 
tural  dynamics  and  control  technology  base  for  LSS.  To  this 
end  we  have  presented  an  assessment  of  both  ACOSS  modeling 
activities  and  control  system  analysis  activities. 

We  have  shown  that  simple  analysis  tools  can  be  utilized  to 
advantage  in  the  LSS  arena.  Through  the  application  of  simple 
and  understandable  tools,  a  "high  order"  structural  model  was 
developed  and  oresented  herein.  Also  "simple"  analysis  tools 
were  presented  for  analysis  of  control  system  stability  and 
performance  characteristics,  both  In  the  continuous  and  sampled 
domains.  The  Insight  and  information  gained  through  application 
of  these  techniques  was  employed  in  participation  with  all  ACOSS 
Contractors  throughout  the  term  of  this  contract. 
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Attn*  <CFf/Capt  dajewski 
F.O.  Box  92960 
rtorldway  Fostal  Center 
Lo.>  Anjelel,  CA  90009 
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Grumman  Aerospace  Oorp 
At».n*  or.  A.  Mendelson 
SoJth  Oyster  day  Poau 
be  t  hpa*  e ,  N7  11/14 


Grumman  Aerospace  Gorp 
Attn*  ..ir.  Art  oertapelle 
Plant  26 

Bechpge,  Hi  11/14 


Jet  Propulsion  Laboratory 
Attn*  .«r.  D.B.  Scnaechter 
48J0  Oak  Grove  Drive 
Pasadena,  GA  91103 


Mli/Lincoln  Laboratory 
Attn*  S.  Wright 
P.O.  box  73 
Lexington,  MA  021 /3 


Mli/Lincoln  Laboratory 
Attn*  Dr.  0.  Hylanu 
P.o.  dox  73 
Lexington,  MA  021/3 


Mif/Lincoln  Laboratory 
Attn*  ur.  A*  Smith 
P.  ).  Box  13 
Lexington,  MA  02173 


Control  Dynamics  Go. 

Attn*  Dr.  Sherman  Seltzer 
22i  Hao t  Slue  Square,  Suite  lb 
Huntsville,  A1  35dJ I 


Lockheed  Space  Missile  Gorp. 

Aten*  A.  Woods,  Blog  130  Organ  62-E6 
3460  Hi  11 view  Ave 
Paio  Alto,  CA  94304 


Lockheed  Missiles  Space  Co. 
Aten*  Mr.  Paul  Williamson 
32b I  Hanover  St. 

Paio  Alto,  CA  94304 
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General  dynamics 
Attn*  Ray  dalstenoerg 


CQivair  division 
SOu  l  Keary  Villa  Rd 
Sail  Diego,  CA  92123 


sri 

AUn*  Mr.  R.C.  Strouu 
20J66  Stevens  Create  bIv.j, 

Cupertiono,  CA  95014 


NASA  Langley  Research  Ccr 
Attn*  Ur.  U.  Horner 

Langley^ 'station  Bldg  I293B  M/s  230 
Hanpton,  VA  23665 

NASA  Johnson  Space  Center 
Attn*  rfooert  HUand 
Ms.  c:A 

riojston ,  IX  7 /Odd 


fflcUonaid  Douglas  Uorp 

At..n*  Mr.  Read  Johnson 

Douglas  Missile  Space  Systems  Div 

53ul  tiulsa  Ave 

Huntington  Beach,  CA  92607 

Integrated  Systems  Inc. 

Attn*  Ur.  Narendra  Gupta 
151  University  Ave. 

Suite  ‘tOJ 

Palo  Alto,  CA  94301 

boeing  Aerospace  Company 
Attn*  Mr.  Leo  Cline 
P.O.  Box  3999 
Seattle  ,  NA  98124 
MS  8  .*-23 

TRw  Dei ense  Space  Sys  Group  Inc. 

Atv.n*  Ralph  Iwens 

Bldg  d2/2054 

One  Space  Park 

ReJondu  Beach,  CA  90278 

TR*. 

Attn*  Mr.  Len  Pincus 
Bldg  H-d,  Room  2031 
Reuondo  Beach,  CA  902/8 


Department  of  the  navy 
Attn*  Ur.  K.T.  Alfrienu 
Naval  Research  Laooratory 
Co-e  /y20 

Washington,  OC  203/6 
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Airesearch  Manur.  Co.  of  Calif. 
Attn*  Mr.  Oscar  Buchmann 
2525  me st  190th  St. 

Torrance,  CA  9U509 


Analytic  Decisions,  Inc. 
Attn*  Mr.  Gary  Glaser 
|4ul  mi  Ison  Blv. 
Arlington,  VA  22209 


Anilyuc  Oecisions.  Inc. 
Attn*  Mr.  Richard  mollicone 
53b6  me st  Rosecrans  Ave 
Sui te  203 

Lawndale,  CA  22209 

Center  for  Analysis 
Attn*  Mr.  jim  Justice 
13  Corporate  Plaza 
Newport  Beach,  CA  92660 


General  Research  Crop. 
Attn*  Mr.  G.  R.  Curry 
P.o.  Box  3587 
Santa  Barbara,  CA  93105 


General  Research  Corp 
Aten*  Mr.  Thomas  Zakrzewski 
7655  Old  Sprinyhouse  Road 
McLean,  VA  22101 


Institute  of  Defense  Analysis 
Aten*  Dr.  rlans  Nolrhard 
40J  Army  Navy  Drive 
Arlington,  VA  22202 


Karman  Sciences  Corp. 

Aten*  Or.  maiter  t.  Rare 
1500  Garden  of  the  Gods  Road 
P.O.  Box  7463 

Colorado  Springs,  CO  80933 

MR j  ,  Inc. 

10100  caton  Place 
Suite  300 
Fairfax,  VA  22030 
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Photon  Research  Associates 
Attn*  mr.  J i*n  Myer 
P.o.  Box  13  Id 
La  Jolla,  CA  92033 


Rockwell  international 

Attn*  Russell  Loftnian  (Jpace  Systems  Oroup) 
(Mail  Code  -  SL56) 

12214  Lakewooa  Ulvd. 

Uowney,  CA  90  24  1 

Science  Applications,  Inc. 

Attn*  Mr.  Richaro  Ryan 
3  P  re st  on  Court 
tied  for-,  4  A  01/  30 


U.J.  Army  Missile  Command 
Attn*  UR5MI-R AS/Mr.  Hrea  Haak 
Redstone  Arsenal,  aL 


Naval  electronic  Systems  Command 
Aten*  Mr.  Charles  Oood 
RMe_IOod4 
National  Center  I 
Washington,  DC  20360 

Naval  Research  Laoora6tory 
Attn*  or.  John  McCallum 
cOfPO 

45-j-j  Overlook  Ave,,  SR 
Waging  ton,  DC  203/5 

U.S .  Army/UARCOM 
Attn*  Mr.  dernie  Chasnov 
A  Me  Bl-g 

50J I  tisenhower  Ave 
Alexan-ria,  VA  22333 


Honeywell  Inc. 

Attn*  br.  Thomas  B.  Cunningham 
Attn*  or.  Michael  F.  Barrett 
2600  Rid  jway  Parkway  MN  I  /-2375 
Minneapolis,  MN  55413 
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MISSION 

of 

Rome  Air  Development  Center 


RAVC  plant  and  execute*  retearch,  development,  teat  and 
*  elected  acgultition  program*  In  support  of  Command,  Control 
Communication*  and  Intelligence  (C3I)  actlvitie*.  Technical 
and  engineering  twpport  within  area*  of  technical  competence 
x*  provided  to  ESP  Program  Office*  (POa)  and  other  ESP 
element*.  The  principal  technical  miction  area*  are 
communication*,  electromagnetic  guidance  and  control,  tur- 
verllance  of  ground  and  aero*paae  object*,  intelligence  data 
collection  and  handling,  Information  tyttem  technology, 
lonotpheric  propagation,  tolid  Atate  tclence*,  microwave 
phytic*  and  electronic  reliability,  maintainability  and 
compatibility. 


